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GTOufPED- 

ABSTRACT 


The  knowledge  of  wave  refraction  is  important  in  many  studies.   The 
rapid  and  relatively  easy  gaining  of  this  knowledge  is  made  possible  by 
the  use  of  the  modern  high-speed  digital  computer.   Large  numbers  of 
spectral  periods  and  incoming  directions  are  easily  investigated,  and 
immediate  results  are  obtained  by  use  of  the  plot  of  the  wave  crest  re- 
fraction from  the  computer.   This  program  presents  the  wave  crest  refrac- 
tion pattern  of  numerous  wave  ray  points  rather  than  the  single  ray 
following  technique.   Its  use  is  valuable  in  amphibious  operation  plan- 
ning, and  in  other  studies  where  the  distribution  of  wave  energy  along 
the  shore  is  desired  for  the  many  periods  of  the  wave  spectrum. 
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1.   Introduction 

The  knowledge  of  the  refraction  pattern  as  an  ocean  wave  approaches 
the  shore  is  necessary  in  the  determination  of  the  energies  of  the  wave 
in  the  surf  zone.   Wave  energy  is  related  in  turn  to  many  near  shore 
processes  or  operations,  such  as  beach  erosion,  sediment  movement,  or 
amphibious  landing  operations. 

While  this  investigation  is  slanted  toward  the  needs  for  amphibious 
operations,  the  knowledge  and  methods  used  are  amenable  to  any  users  who 
desire  the  pattern  of  wave  crest  or  refraction,  or  an  estimate  of  the 
ratio  of  the  wave  height  at  a  point  to  the  deep-water  wave  height  (H/H  ) 
for  other  purposes. 

Early  work  on  the  subject  was  concerned  mostly  with  hand  calculation 
of  the  various  determining  parameters  and  graphical  method  for  construc- 
ting the  refraction  diagrams  |_lj  ^2j  .   Later,  in  early  1962,  attempts  were 
made  to  evaluate  wave  refraction  using  high-speed  digital  computers  [3j. 
The  most  usable  results  published  were  those  obtained  by  Harrison  and 
Wilson  [k\ .      While  the  work  of  Griswold  and  Nagle  as  well  as  that  of  Har- 
rison and  Wilson  gives  worthwhile  results,  both  studies  are  very  limited 
with  respect  to  the  ratio  of  work  input  to  the  results  obtained.   They 
are  concerned  with  following  one  ray  at  a  time  from  deep  water  to  the 
shoreline.   To  get  a  complete  refraction  diagram  many  inputs  to  the  com- 
puter are  required,  and  hand  plotting  of  the  results  are  necessary. 

The  procedure  outlined  in  this  paper  follows  a  wave  crest  composed 
of  a  number  of  ray  points  to  the  shoreline  with  immediately  available 
results.   No  hand  plotting  or  further  program  input  is  required  for  the 
rays  with  one  period  and  the  same  initial  deep-water  direction.   By 
changing  only  the  wave  period  or  wave  direction,  the  spectrum  of  periods 


and  range  of  angles  can  be  investigated  rapidly  with  little  additional 
time  involved. 

All  of  the  computer  programs  for  the  investigation  of  wave  refraction 
utilize  the  same  basic  procedure:   following  the  wave  ray  (orthogonal)  to 
the  shoreline.   The  new  program  uses  a  depth  field  as  an  interpolating 
surface  rather  than  the  velocity  field  of  [3J  and  [4J .   Several  interpo- 
lation surfaces  have  been  investigated  to  represent  the  velocity  or  depth 
field  used  [4j .   This  program  utilizes  a  quadric  surface  for  interpola- 
tion of  the  depth  field. 

The  present  study  uses  an  input  of  deep-water  wave  period  and  wave 
direction  to  a  computer  program  for  determining  wave  refraction  in  a 
method  such  as  proposed  by  Munk  and  Arthur,  which  employs  the  wave  para- 
meters listed  in  Appendix  IV  (_5 J  .   An  additional  subroutine  of  the  pro- 
gram computes  the  coefficient  of  refraction,  K,  and  the  ratio  of  wave 
heights,  H/H  .   The  latter  two  values  are  recorded  along  with  the  X  and 
Y  coordinates  of  each  point  along  the  wave  crest.   As  a  rapid  means  of 
viewing  the  refraction  pattern  of  the  wave  crest,  a  graphical  output  is 
included  which  contours  every  third  wave  crest  computed.   Other  parameters 
may  be  included  in  the  output  at  the  user's  discretion. 


2.  Method. 

The   first   step   in  the  utilization  of  the  program   is   the  construction 
of  a  grid  of  depth  values   for  the  desired  area.      This   grid  must   include 
starting  points   in  deep  water  for  all  rays   to   be   followed,    such  that   the 
ratio  of  the  depth  to   the  deep-water  wave   length,   d/L0,    is  greater  than 
0.5.      Since  the  program  is   arranged  to   follow  the  waves   from  deep  water 
to  the   shore,    the  grid  origin  must   be   in  deep  water.      The   convention   is 
that  the  X-axis  will   be  positive   and   increasing  toward  the   shore  while 
the  positive  Y-axis    is    90   degrees   to   the   left  of   the  X-axis    as   shown   in 
Figure   1.      The  grid   interval    is    selected   such  that,    in  a  given  cell,   the 
bottom  contours   are  reasonably  parallel  to  one   another. 

Actual  or   interpolated  depths   at  the  grid   intersections   are  recorded 
to   the  best   accuracy  available   from  the  chart,    and   all   actual   depths   are 
made  positive.      Extrapolated  depth  values   are  continued  on   land   for  two 
grid   units   from  the   shoreline   and   are  made   negative.      For  depths  on  the 
shoreline   itself,    zero   is   used.      Depths  on   land  outside  of  two   grid  units 
from  the   shoreline  are  made   some  arbitrary,    negative,    non-zero  value. 
The  procedure   for   assigning   negative  values   for  the   nearshore   land  depths 
is   required   in  the  fitting  of   a  surface  to   the   localized  depth  values. 
The   zero   land  depths   are  used   for  contouring  the   shoreline  on  the  output 
graph. 

3.  Example  of   Input. 

As   an  illustration  of  the   procedure  and  to   test  the  results,    the 
southern  portion  of  Monterey   Bay,    California,   was   used.      The  depth  grid 
was   selected   so   that  the  origin  of  the   18  rays  would   be   positioned   in 
water  depths  greater  than   1,024  feet  where  d/L     =0.5   for  a  period  of   20 
seconds.      The  direction  toward  which  the  rays   proceed    is    125°  TRUE,   making 


an  angle  with  the  X-axis  (oriented  at  90°  TRUE)  of  -  35°.  The  grid  inter- 
val selected  was  1,500  feet,  and  the  depth  input  was  in  fathoms  to  facil- 
itate the  chart  reading.   The  depth  values  were  determined  from  the  U.S. 
Coast  and  Geodetic  Survey  chart  540  3. 

Other  parameters  that  were  required  for  input  were  the  X  and  Y  values 
of  the  starting  point  for  the  first  ray.   The  values  of  the  ray  separa- 
tion parameter,  g  ,  are  input  as  Bl  and  B2,  and  are  always  equal  to  1.0 
in  deep  water.   The  period  of  the  waves  and  the  angle  with  respect  to  the 
X-axis  are  entered  as  T  and  Al ,  respectively.   The  following  parameters 
are:  NOR,  the  number  of  rays  that  are  being  followed;  DIST,  the  distance 
between  these  rays;  TIME,  the  time  interval  for  advancement  of  the  wave 
front;  GRID,  the  grid  interval;  MM,  the  number  of  grid  points  in  the  X 
direction;  and  NN,  the  number  of  grid  points  in  the  Y  direction.   The 
angle  of  the  wave  direction  is  in  relation  to  the  X-axis  and  is  the  di- 
rection toward  which  the  waves  are  moving.   It  can  be  positive  or  negative 
with  respect  to  the  X-axis. 
4.   Computer  Operations. 

The  computer  first  reads  the  parameter  values,  then  the  depth  grid 
as  a  column  of  Y-values  for  a  constant  X-value.   The  depths  are  converted 
to  feet  immediately.   The  subroutine  DEPTFUN  is  called  to  compute  the 
depth  at  the  first  point  by  fitting  the  closest  nine  grid  point  depths 
around  the  starting  value  to  a  quadric  surface  by  the  least  squares  meth- 
od, using  an  equation  of  the  form: 

DEPTH    =  Af  +  AZX  +  A3r  +  A+XY  +  A5X2  +  AL  Y* 

where  DEPTH  is  the  value  of  the  depth  at  that  point,  A..  „    are  constant 

J. ,2.  . n 

coefficients,  and  X  and  Y  are  the  distance  values  for  the  point.   Each 


time  a  new  depth  is  encountered  the  surface  of  best  fit  is  computed  from 
the  surrounding  nine  grid  points.   Also  included  in  the  DEPTFUN  subrou- 
tine is  the  computation  of  the  wave  velocity  at  that  point  for  the  given 
depth.   This  is  done  by  an  iteration  process  using  the  common  result  from 
solving  the  wave  equation  as  used  in  H.O.  Pub.  234  [2]  : 


T 


-"   "JJ 


Zr  (  TCJ 

where  C  =  wave  velocity,  g  =  acceleration  due  to  gravity,  T  =  wave  per- 
iod, and  d  =  water  depth.   As  with  other  wave  refraction  investigations 

1 

it  is  assumed  that  the  wave  velocity  is  a  function  only  of  water  depth 
and  wave  period.   Other  factors  such  as  bottom  friction,  percolation, 
currents,  reflection,  and  wind  are  considered  as  not  affecting  the  re- 
fracting waves. 

The  wave  ray  is  moved  to  the  next  point  by  solving  for  the  ray  curva- 
ture and  projecting  the  ray  forward  in  the  time  interval  specified  at  the 
speed  calculated  for  the  point. 

The  curvature  is  calculated  using  the  expression  [5J: 


i  faA$    -     c~4$] 


where  FK  =  ray  curvature,  and  A  =  approach  angle. 

To  determine  the  values  of  X   ,  ,  Y  . ,  A„  , ,  ,  and  FK  ,  for  the  suc- 

n+17  n+17   n+17       n+1 

ceeding  point,  an  iteration  process  is  used  to  solve  the  equations  [3J: 

A  A     =  (^    +   f^()  1>/Z 
4h+,    =    An   +  AA 

A       =  (A*  +  A^)/z 
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/  r        /     f     D     C6S>    A 

Yh+,       =     Yn    +  t>    sin  A 


where  D  =  the   incremented  distance    (D  =  CT)   between  points   n  and  n+1. 
At   the   point   n+1,   the  value  of   Beta   is   calculated  by   solving  the 
second-order,   non-linear  differential   equation    \5 J : 


$& 


b% 


2 


+   P'ff  +  1$    =  ° 


P*    ~    cc,  A    £*     .    ,„A    i& 


where 

The  above  equation  is  solved  by  the  finite  difference  method  [3  J.   This 
results  in  an  equation  for  the  Beta  value  at  the  n+1  point  in  terms  of 
the  Beta  values  at  the  two  previous  points.  The  equation  to  be  solved 
is  then: 

Z  +  pl> 
where  p,  q  and  D  are  defined  above,  and  u      and  p„  are  the  Beta  values 
of  the  two  previous  points.   It  is  to  be  noted  that  the  calculation  of 
Beta  is  made  for  the  n+1  point.   In  the  BETA  subroutine  the  coefficient 
of  refraction  is  calculated  by  the  relation: 


K~{T7p 


where  K  =  the  refraction  coefficient. 

To  determine  the  ratio  of  wave  height  at  the  point  to  the  deep- 
water  wave  height,  the  shoaling  factor  has  been  approximated  for  values 
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of  C/C0  <  1  by  a  curve   such  that    [2]  : 

//,  =  3.2S/9  -  iz.ai&fy  +  Z8.6ff)k(£)2  -  a?,  ?zs7  (0 

4 

+      ll.68/f(j;) 

La 

As  C-.,   the  deep-water  wave  velocity,    is   a  function  of   period   alone 
and  C    is   known  for  any  point    (X,Y),   the  equation  may  be  evaluated   at   any 
or  every  point  at  which  the  calculations   are  made  for  C.      The  ratio  of 
wave  heights    (H/FL)    is   found   from  the  equation  [2j: 

5.  Output   from  Computer. 

The  printed  output  from  the  computer  consists  of  the  X  and  Y  values 
of  each  position  after  advancement  by  the  increment  of  time  specified. 
This  is  in  units  of  yards  for  ease  in  hand  plotting  on  the  charts.  The 
coefficient  of  refraction  and  the  wave  height  ratio  for  each  point  are 
given. 

Included  in  the  program  output  is  a  graphical  plot  of  the  wave  crests, 
which  is  programmed  for  the  Calcomp  160  computer  system,  utilizing  the 
DRAW  subroutine  of  Appendix  II.   The  first  and  every  third  crest  there- 
after are  plotted  and  points  connected  depicting  the  wave  crest.   The 
scale  is  such  that  a  true  representation  of  the  wave  crest  is  presented. 
From  this  graph  and  if  desired,  hand  plotting,  areas  of  convergence  and 
divergence  are  easily  seen.   By  knowing  the  number  of  the  crest,  the  para- 
meters of  refraction  coefficient  and  wave  height  ratio  can  be  found  from 
the  printed  data. 

6.  Program  Development. 

To  calculate  the  coefficient  of  refraction  and  the  wave  height  ratio 
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it  is  necessary  that  the  interpolation  surface  for  the  depth  values  be  of 
the  second  order,  at  least,  so  that  the  second  partial  derivatives  would 
be  available  for  the  computation  of  these  values.   Consequently,  a  quadric 
surface  was  used  for  the  representation  of  the  depth  values  at  the  closest 
nine  grid  points  to  the  position  at  which  the  values  were  desired.   The 
partial  derivatives  of  the  surface  in  the  X  and  Y  directions  were  used  to 
evaluate  the  partial  derivatives  of  the  velocity  function  as  proposed  by 
Harrison  and  Wilson  |4| .   However,  a  power  series  representation  of  the 
hyperbolic  tangent  was  used  to  evaluate  the  velocity  derivatives  rather 
than  the  method  used  by  Harrison  and  Wilson,  as  shown  in  Appendix  III. 

The  calculation  of  the  Beta  parameter  is  required  for  determining 
the  refraction  coefficient.   The  finite  difference  method  is  used  in  this 
program.  However,  the  use  of  this  method  requires  that  the  distances  be- 
tween the  points  at  which  the  equation  is  being  solved  be  equal.   In  the 
method  presented  here  where  the  distance  is  a  function  of  the  velocity  at 
each  point,  this  does  not  hold  exactly  true.  The  difference  in  the  veloci- 
ties at  the  successive  points  is  on  the  order  of  one  foot  per  second  due 
to  the  shallow  contour  gradient  in  the  particular  area  of  interest.  This, 
of  course,  will  not  hold  true  for  all  cases.  A  better  method  would  be  to 
solve  the  second-order  equation  in  Beta  by  the  Kelvin  method  which  re- 
quires only  the  distance  between  the  point  n  and  n+1,  which  is  readily 
obtainable. 
7.   Discussion  of  Results. 

The  printed  results  are  shown  in  Appendix  I  for  several  of  the  wave 
crests.   The  X  and  Y  values  from  this  type  of  presentation  were  plotted 
by  hand  as  in  Figure  1.   This  plot  shows  areas  of  marked  divergence  along 
the  shore.   When  this  figure  is  compared  with  Figure  2,  the  result  of  the 
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graphical  method  of  H.O.  Pub.  234,  little  difference  is  noted  suggesting 
that  both  the  computer  method  and  the  graphical  method  produce  similar 
results.   This  is  the  aim  of  the  investigation.   Other  comparisons  were 
made  using  different  directions  and  periods  with  comparable  results. 

In  Figure  1,  the  seventh  and  eighth  wave  rays  are  seen  to  cross  and 
continue  to  the  shore.   This  crossing  is  attributed  to  the  bottom  contours 
of  ray  eight  when  the  ray  first  reaches  shallow  water  (d/Ln  =  0.5).   There 
is  a  small  but  steep-sided  indentation  in  the  otherwise  gradual  slope  of 
the  area.  This  causes  the  ray  to  bend  due  to  the  steep  depth  gradient 
which  results  in  convergence  with  ray  seven  below  it.   From  this  point 
to  the  shore  the  refraction  is  similar  to  that  of  ray  six.   Figure  3,  the 
computer  drawn  plot,  shows  all  of  the  wave  crests  from  the  first  advance- 
ment to  the  shore.   The  ray  crossing  is  evident. 

The  values  of  the  coefficient  of  refraction  and  the  wave  height  ratio, 
Appendix  I,  show  little  refraction  for  crests  one  and  two  as  expected, 
since  the  waves  for  the  most  part  are  in  deep  water.   However,  later  values 
of  these  parameters  do  show  the  effects  of  the  refraction  seen  in  Figures 
1  and  3.   The  values  of  these  parameters  require  verification.   The  values 
do  appear  to  be  qualitatively  reasonable  for  the  refraction  encountered, 
and  when  compared  to  the  values  estimated  by  the  technique  of  measuring 
the  distance  between  orthogonals. 

It  is  estimated  that  the  time  necessary  to  construct  the  grid  and 
record  the  depth  values  is  three  to  four  hours  depending  on  the  size  of 
the  area  desired  and  the  gradient  of  the  bottom  contours  in  that  area. 
A  shallower  gradient  requires  a  larger  area  to  ensure  that  the  wave  rays 
start  in  deep  water.   The  time  required  to  punch  the  data  cards  and  to 
check  the  results  is  a  function  of  the  experience  of  the  operator.  The 
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computer  time  is  approximately  five  minutes  for  the  computation  of  the  re- 
sults in  this  paper,  Appendix  I,  plus  an  additional  five  minutes  to  draw 
the  graph.   This  time  can  be  reduced  drastically  if  the  program  as  com- 
piled by  the  computer  is  put  on  tape  and  input  in  that  form.   The  compiler 
takes  over  half  of  the  computer  time  at  present.   The  time  is  spent  con- 
structing the  grid  of  depth  values  only  once.   From  this  grid  all  of  the 
various  wave  periods  and  directions  can  be  investigated  without  further 
effort  of  constructing  another  grid  as  is  necessary  in  the  hand  plotting 
of  the  refraction  diagram. 

8.  Conclusion. 

As  discussed  in  the  preceeding  section,  the  results  of  the  refraction 
diagram  agree  well  with  that  of  the  hand  drawn  method,  so  that  the  plots 
received  from  the  computer  can  be  used  with  as  much  confidence  as  those 
drawn  manually.   The  decided  advantage  is  that  the  computer  product  can 
be  obtained  many  times  faster. 

The  coefficient  of  refraction  and  the  wave  height  ratio  as  noted 
appear  qualitatively  correct  and  can  be  used  with  the  reservation  that 
the  values  have  not  yet  been  verified. 

9.  Recommendations. 

Further  development  of  programs  similar  to  this  will  require  a  bet- 
ter representation  of  the  bottom  contours  for  more  accurate  results.   The 
shallow-water  depth  values  need  to  be  very  accurate  for  the  proper  deter- 
mination of  all  the  parameters,  since  the  depth  determines  the  velocity 
of  the  wave  at  each  point,  and  the  other  parameters  are  functions  of  the 
velocity. 

An  accurate  method  within  the  approximations  made  is  required  for 
the  calculation  of  the  Beta  parameter.   The  analytic  approach  of  Kelvin's 
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method  appears  to  be  the  more  accurate  method  from  the  consideration  of 
the  necessary  assumptions  involved  with  the  finite  difference  method.   In 
the  case  involving  steep  gradients  of  the  bottom  contours,  the  finite  dif- 
ference method  as  used  here  may  prove  inadequate.   The  time  factor  in  com- 
pleting this  paper  prevented  the  successful  completion  of  programming  the 
Kelvin  method  as  desired. 

Testing  of  the  results  is  required  to  substantiate  the  procedure  as 
a  prediction  method.  Involved  aerial  photography  is  required  to  follow 
the  wave  crest  refraction  to  the  shoreline  where  the  accurate  prediction 
is  most  desired.  Wave  gages  and  other  devices  may  be  used  to  obtain  the 
necessary  data  for  the  verification  of  the  wave  height  values  as  predic- 
ted by  the  computer.  A  very  intensive  project  will  be  required  in  order 
to  verify  the  predicted  values. 
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X-SCft_E  -   1.QQE+8H  UHITS/INCH 
Y  -SCALE  -  L08E+»i  UMITS/JMCH. 

STOUPPE  WflUE  REFRACTION  PROGRAM 

ANGLE  =    -35  DEG,   PERIOD  =  26  SEC. 

Figure  3.   Computer  drawn  plot  of  wave  crests 
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APPENDIX  I 

Computer  Program  for  Wave  Refraction  using  Fortran  63 

on  the  C.D.C.  1604  Computer 

Program  Title;   WAVREFR 

Input  Variables: 

XV,  YV  X  and  Y  values  of  the  first  wave  ray  to  be  computed  (feet). 

Bl,  B2 Values  of  the  Beta  coefficient  (always  equal  to  one  in  deep 

water) . 

T  Wave  period  (seconds). 

Al  Angle  measured  from  the  direction  of  increasing  X  along  the 

X-axis  in  the  direction  of  travel  of  the  wave  crest  (degrees). 

NOR Number  of  rays  that  will  be  followed. 

DIST  Perpendicular  distance  between  rays  in  deep  water  (feet). 

TIME  Time  interval  used  for  advancing  the  wave  front  (seconds). 

GRID  Distance  between  grid  points  of  depth  values  (feet) . 

MM,  NN  Number  of  grid  points  in  the  X  and  Y  directions,  respectively. 

DEP  (I, J)...  Depth  values  in  the  grid  (fathoms). 

Output  Variables; 

X,  Y  X  and  Y  coordinates  for  a  given  wave  crest  and  ray  number 

(yards) . 

COREFR   Coefficient  of   refraction  for  the  wave  at   the   point  X,Y. 

HHO Ratio  of  wave  height   to   the  deep-water  wave  height   at   the 

point    X,Y. 

NGO  Indicates  that  the  ray  has  terminated  (=  1),  or  is  continuing 

(=  2). 

Variables    in  Common;    (not   previously  defined) 

FK,   FFK   ....    Values   of   ray  curvature. 

CXY,  CXX   . . .    Values   of  wave  velocity. 

PAR Constant    =   gT/2ir. 
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BAR Constant  =  2  TT  /T. 

MAX  Number  of  wave  crest  (=  1  for  first  crest). 

DEPTH  Calculated  depth  at  point  X,  Y. 

NORA Number  of  ray  being  followed. 

ID Indicates  ratio  of  d/L0- 

=  1  for  d/Lg  >  0.5; 

=  2  for  d/L0  £  0.5. 

MIT  Designates  whether  the  last  two  curvature  estimates  for  a 

given  X,  Y  are  less  than  .0001,  whether  the  13th  and  15th 
estimates  are  less  than  .0001,  or  whether  neither  of  the 
above  is  true  (MIT  =  1,  2,  3,  for  the  respective  cases). 

NIP  Used  to  determine  the  number  of  wave  crest  to  be  graphed. 

XXX,  FX  ....  Values  of  X  used  in  the  graphing  and  printout. 

YYY,  FY  Similar  to  uses  for  XXX  and  FX. 

NNGO  The  number  of  rays  that  have  stopped. 

XX,  YY,  AA,  A 

Intermediate  values  of  X,  Y,  and  A. 

Summary  of  Program: 

The  program  reads  the  input  variables,  then  the  depth  grid.   PAR  and 
BAR  are  computed  for  the  wave  period;  A  is  converted  to  radians;  and  MAX 
and  NORA  are  set  equal  to  one.   An  initial  value  of  CXY  is  found  by  test- 
ing the  wave  period.   Control  is  transferred  to  the  RAYN  subroutine.   When 
RAYN  has  determined  the  parameters  at  the  first  and  second  points  along 
the  ray,  NORA  is  increased  one,  the  XX  and  YY  values  computed  for  the  next 
ray,  and  then  RAYN  is  called  again  to  compute  the  parameters.   This  pro- 
cedure is  followed  until  NORA  equals  NOR,  the  number  of  rays.   The  graph 
is  drawn  for  the  first  wave  crest  and  the  values  printed.   MAX  is  in- 
creased by  one,  and  the  procedure  started  over  until  all  the  rays  have 
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stopped  by  either  going  off  the  grid,  hitting  the  shoreline,  or  not  having 
the  values  of  the  curvature  converge.   When  all  the  rays  have  stopped  NNGO 
is  equal  to  NOR.   The  final  graph  is  drawn  which  is  the  contour  of  the 
area  utilizing  the  zero  values  of  the  depth  grid  and  the  program  is  com- 
plete. 
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Subroutine  Title;   RAYN 
Variables  o£  Subroutine; 

COREFA  Intermediate  value  of  COREFR. 

HHA Intermediate  value  of  HHO. 

FKBAR Curvature  used  to  obtain  DEL  A. 

NOGO Storage  value  of  NGO. 

AAA Storage  value  of  A. 

Summary  of  Subroutine; 

RAYN  calls  DEPTFUN  to  obtain  CXY,  PDPX,  PDPY,  PDDPXY,  PDDPXX,  and 
PDDPYY.   NGO  is  set  equal  to  two  if  DEPTH  is  not  zero,  otherwise  NGO  is 
set  equal  to  one.   RAYN  tests  NGO  of  the  ray  for  the  MAX  calculation  to 
determine  whether  the  ray  will  continue  or  not.   If  NGO  =  2,  KFUNCT  is 
called  to  obtain  the  value  of  the  curvature,  FK.   If  NGO  =  1,  control  is 
returned  to  WAVREFR.   If  the  ray  is  continuing  RAYN  calls  MOVE  to  project 
the  ray  to  its  next  position.   If  the  ray  is  not  stopped  in  the  MOVE  sub- 
routine, RAYN  next  calls  BETA  to  calculate  the  values  of  COREFR  and  HHO. 
As  a  final  step  RAYN  stores  the  values  of  FK,  COREFR,  HHO,  NGO,  X,  Y, 
PDPC,  and  AAA  for  use  when  the  ray  is  again  projected. 
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Subroutine  Title:   DEPTFUN 

Variables  of  Subroutine; 

KER Indicates  errors  in  the  solution  of  the  simultaneous  equa- 
tions for  the  quadric  surface: 

=1   no  errors  indicated  in  the  solution; 

=  2   indicates  that  the  matrix  of  values  being  solved  is  singular 
or  nearly  singular. 

ALO  Deep-water  wave  length. 

DLO Ratio  of  DEPTH  to  ALO. 

PDPX,  PDPY,  PDDPXX,  PDDPXY,  PDDPYY 

First  and  second  partial  derivatives  of  the  depth  at  X,  Y 

with  respect  to  the  X  and  Y  directions- 

Summary  of  Subroutine: 

The  subroutine  first  determines  the  values  of  the  closest  grid  point 

by  truncating  the  values  of  XX/GRID  +1.5  and  YY/GRID  +  1.5  to  give  the 

correct  value  of  the  grid  point.   The  quadric  surface  is  fitted  to  the 

nine  values  of  depth  surrounding  this  calculated  point,  using  the  least 

squares  method.   DEPTH  is  found  by  evaluating  the  quadric  equation  at  X 

and  Y  values  of  the  point.   If  DEPTH  is  positive,  the  wave  velocity,  CXY, 

is  solved  for  by  an  iteration  procedure  of  the  equation  described  in  the 

test  using  the  principles  of  Appendix  IV.   The  various  partial  derivatives 

are  computed  by  evaluating  the  equation  described  in  Appendix  III.   The 

method  for  the  solution  and  the  program  for  solving  the  six  simultaneous 

equations  resulting  from  the  least  squares  method  was  written  by  C.  B. 

Bailey  and  Mary  Haynes  of  the  USNPGS  computer  center  programming  staff. 
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Subroutine  Title:      KFUNCT 


Variables  of   Subroutine: 


PDPC,   PDDPCC 


First  and  second  partial  derivatives  of  the  depth  with 

respect  to  the  wave  velocity. 

PCPX,  PCPY  ..  First  derivatives  of  the  wave  velocity  with  respect  to  the 
X  and  Y  directions  calculated  from  relations  in  Appendix 
III. 

FK Curvature,  of  the  ray  at  the  point  X,Y. 

Summary  of  Subroutine: 

The  subroutine  calculates  the  values  of  PCPX  and  PCPY,  and  the  curva- 
ture at  the  point  X,  Y  from  the  equations  of  Appendix  III  and  the  text, 
respectively.   The  curvature  is  calculated  only  if  LO  =  2 ,  where  the  wave 
is  in  shallow  water. 
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Subroutine  Title;   MOVE 
Variables  of  Subroutine; 

FFKK  Storage  value  of  FKBAR. 

DEL  D  Increment  of  distance  advanced  (DEL  D  =  T  CXY). 

DEL  X XX  -  X. 

DEL  Y YY  -  Y. 

DEL  A  AA  -  A. 

ABAR (A  +  AA)/2. 

Summary  of  Subroutine; 

MOVE  determines  the  X  and  Y  values  of  the  next  point  on  the  wave  ray 
by  an  iteration  process  involving  the  curvature,  the  incremented  distance, 
and  the  angle  of  the  ray.   The  iteration  is  continued  until  the  curvature 
estimates  vary  no  more  than  .0001  from  one  another.   Then  the  values  of 
XX,  YY,  AA,  and  FKK  are  accepted  for  the  new  point.   If  the  difference  is 
greater,  FKKP  is  set  equal  to  FKBAR  and  the  current  FKBAR  is  used  to  obtain 
another  set  of  values.   If  the  iteration  process  stops  before  15  itera- 
tions, MIT  =  1.   If  the  cycle  stops  at  15  iterations,  and  the  difference 
between  FKBAR  and  FKKP  is  less  than  .0001,  MIT  =  2,  and  FKBAR  is  defined 
as  (FKBAR  +  FKKP)/2  for  the  last  determination  of  XX,  YY,  AA,  and  FKK 
values.   If  the  difference  is  greater  than  .0001  after  15  iterations, 
MIT  =  3,  and  the  ray  is  stopped.   Control  is  sent  back  to  RAYN.   Inside 
the  iteration  loop  DEPTFUN  determines  the  parameters  of  CXY,  PDPX,  and 
the  other  derivatives,  and  tests  the  values  of  X  and  Y  to  determine  if 
they  are  still  on  the  grid. 
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Subroutine  Title;     BETA 

Variables  of   Subroutine; 

PCCPXX,    PCCPXY,    PCCPYY 

Second  partial  derivatives  of  the  wave  velocity  with  respect 

to  the  X  and  Y  directions. 

HSHOL   Shoaling  factor. 

CCO   Ratio  of  the  wave  velocity  to  the  deep-water  wave  velocity. 

DD    Increment  of  distance  along  the  wave  ray  between  points  n 

and  n+1. 

B(i,j)  Values  of  the  Beta  parameter  at  the  three  points. 

Summary  of  Subroutine; 

BETA  calculates  the  values  of  the  coefficient  of  refraction  and  the 
wave  height  ratio  at  each  point  along  the  ray.   The  equations  are  des- 
cribed in  the  text.   If  LO  =  1,  BETA  sets  COREFR  and  HHO  equal  to  one, 
and  the  Beta  value  equal  to  the  previous  Beta  value  on  the  ray.   The 
shoaling  factor  is  calculated  from  an  equation  determined  from  a  poly- 
nomial fit  to  a  curve  of  C/Cn  <  1  as  described  in  the  text. 
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PROGRAM  WAVREFR 
C      WAVE  REFRACTION  PROGRAM  FOR  THE  COMPUTATION  OF  CREST  AND  RAY 
C      REFRACTION  IN  SHAOLING  WATER  BY  LT.  D.  STOUPPF.  OCTOBER  1966. 
C      THE  PROGRAM  IS  WRITTEN  IN  FORTRAN  63  FOR  THE  CDC  1604  COMPUTER 
C      SYSTEM,  INPUT  OF  THE  BINARY  IS  REQUIRED  FOR  THE  DRAW  SUBROUTINF 
C      ON  THE  PRESENT  SYSTEM  AT  THE  US  NAVAL  POSTGRADUATE  SCHOOL. 
C      NOTE  TO  USERS   ALL  INTEGER  INPljT  ON  THE  DATA  CARDS  MUST  BF 

RIGHT  ADJUSTED  FOR  CORRECT   RESULTS. 

USERS  KAY  ENTER  THE  DESIRED  NAMES  FOR  THE  GRAPH  TITLES  IN  THE 

ITITLE  STATEMENTS  AS  DESCRIBED  IN  THF  WRITFUP  ON  THE  DRAW 

SUBROUTINE. 

/DIMENSION  ITITLEC 12)  .XXX ( 100) ,YYY( 100)»FX(100),FY(100) 

COMMON/BLK1/X1100) »Y( 100) ,CORPFR(100) ♦HH0( 100 ) .B ( 3 » 100 ) ,FFK ( 100 ) , 
INOGOUOOWNGO.AAAi  100) 

C0MMGN/eLK2/TiA.CXY.PAR,BAR,TlME.GRID.FK»MAX,N0R,MM,NN,DEPTH,HHA, 
lC0REFA»N0RA,L0tDIST 

COMMON/ BLK3/DEP( 100.10  0) .PDPX .PDPY .PCPX . PCPY .PDDPXX ♦ PDDPYY 
•  WPDDPXY.POOPCC.  POP  C 

READ2.XV.YV.B1.B2.T.A1.NOR.DIST.TIMF.GRID.MM.NN 
2  FORMAT  (2F10.1.2F3.0,F4.0,F4.0,I4,F10.2,F5.1 .F10.1.2I4) 

PRINT  19 

19  FORMAT  {  1H1»5X.2HXV,7X.2HYV,4X,?HB1  » 1X.2HB2  1 3X.»  1HT  »2X  »2HA1  •  ?X  » 
13HNOR.3X.4HDIST.4X.4HTIME.3X»4HGRID.3X,2HMM.2X»2HN'M/) 

PRTttT  2. XV  .YV.B1.B2.T.  A  l.NOR.D  IST.TIME.GR  ID.  MM.  NN 
READ  l.'(  (DEP(  I.J) .J=1.NN) .I=1,MM) 
1  FORMAT  (14F5.0) 

OOZZ    I«1.WM  $  D022  J=1.NN 

22  DEP(I»J)*DEP( l.J)*6. 
MAX*1  $  A1«A1*. 01745329  $  PAR=32.2*T/6 . 28318 5 

LABEL-4H        $  DO  18  M=l,12 
18  1TITLE(M)»8H  $  ITITLE(1)=8H  STOUPPc  S 

ITITLE (  5)  «8>fRACT TON    $  IT  ITL^  (  6  )  =8HPROGRAM 
ITITLE(8)*8HANGLF  =   $  IT  I TLF ( 10 ) =8H  PFPIOD 
ITITLE(9)=8H-35  DEG,    $  I T ITLF ( 11 ) =8H=  20  SFC 
e(l.l)«Bl  $  8<2,1)=B2  $  XX=XV  $  YY=YV  t  NORA=l 
If (T-J0.J15.16.17 

15  CXY«  30.0  $  GO  TO  3 

16  CXY»  50.0  $  GO  TO  3 

17  CXY«80.0  $  GO  TO  3 
5  B(1.N0RA)=B1  $  B(2.NORA)=B2  J  A=A1 

N   •  XV*XV-DIST*SINF(A)  $  YV=  YV+Dl ST*COSF ( A ) $XX=XV 

23  NORA=l 

24  XX*X(NORA)  $  YY=Y(NORA)  $  A=AAA(NORA) 
5  CALL  RAYN(XX.YY.NOR) 

CONTINUE  $  NORA  =  NORA  +  l  $  I F ( Noda-NOR ) 10  ,  10  ,4 
10  !F(MAX-1)5.5.24 
4  IF<MAX-1)20,20.25 

25  IF(NIP-4)21.26,26 

20  CALL  DRAW(N0R. X.Y.I. 2. LABEL. ITITLE. 10000. .10000. ,0,0. 2. 2.9, 9  .1 
1LAST)  *  GO  TO  27 

26  CALL  DRAW(NOR.X.Y.2.2.LABEL,ITTTLE,100O0.,10O00.,O,O,2,2,9,9  ,1 

1LAST) 


S  BAR=6.283185/T 


ITITLE(4)=8HW*Vfr  R^ 


S  A  =  A1 


5  YY=YV  J  GO  TO  3 


WAVFOOOO 
WAVEOOIO 
WAVEO020 
WAVE0030 
WAVE0040 
WAVEO050 
WAVE0060 
WAVE0070 
WAVEOO80 
WAVFOH90 
WAVFOIOO 
WAVE0110 
WAVF0120 
WAVF0130 
WAVE0140 
WAVE0150 
WAVE0160 
WAVE0170 
WAVF0180 
WAVF0190 
WAVF0200 
WAVFO^IO 
WAVF0220 
WAVFA230 
WAVF0  240 
WAVE^250 
WAVE026O 
WAVE0270 
WAVF0280 
WAVEH290 
WAVF03OO 
FWAVr^10 
WAVCn^20 
WAVF0^30 
WAVF0^40 
WAVF03SO 
WAVFrt360 
WAVF037O 
WAVE0380 
WAVFO390. 
WAVE0400 
WAVEO410 
WAVF0420 
WAvEn430 
WAVr^440 
'.fAVc045O 

WAVC^46A 
WAVF0  47O 
WAVF04PO 

tfAVEO*^ 

WAVP0  5HO 
WAVE^MO 
V AVE ^5  20 
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27  CALL  DRAW(N0R.»X»Y»2,  0, LABEL,  I  T  !  TLC  ,  nOO  0  .  ,10000.  i0>0,2»2,9,9  ,  v,   WA.VF.0  53n 
1  LAST)  vAVCr>K4n 

MIP  =  1  .WAVE  0*50 

21  NIP  =  NIP  +  1  S  MAX=MAX+1  S  MIN=MAX-1  WAVE0560 

33  00  32  M=1»N0R  X  FX  (  M  )  =X  (  M  )  /3  •  5  CY (  V1  =Y  <  V  )  /3  •  .      WAVErt*70 
32  CONTINUE      v  WAVEn«>SO. 

34  PRINT  28»*4IN  "'  WAVE0-59O 

28  FORMAT  (//17HNUM3ER  OF  CR-ST  =  14/)  V.'AVpn6*n 
PRINT  6  WAVF0610 

6  FORMAT  (  /4X»1HX»  12X  » 1HY  >8X  »  6HCO&EFR  »3X»'3HHHO  »  7X  >3HNGO,  10X  »  1HX*12X»  '''AVF0  6  2R 
11HY,8X,6HC0REFP,3X  ,  3  mHHG  ,  7X  ,  3^00/  )  WAVE0630 

DO  7  J=l»NOR»2  WAVE'n64* 

PRINT  8,FX(  J)  »F'Y.(  J)  »COREFR(J)  ,mH'Q(  J)  ♦  NOGO  (  J  )  ,  FX  (  J  +  l  )  ,  <^Y  (  J  +  l  )  ,  WAVE.  065* 

1C0REFR(  J+1)  ,HH0  (  J+l  )  ,NOGO<  J  +  l  )  WAVEn6£n 

8  FORMAT  (  2(F10.2,3X,F10.2,'5X,fr5.'',^y,crC#c,3y,I5,7X)/)  VAV€n6?" 

7  CONTINUE  WAWF.068* 
DO  9  I  =  l»NOR  S  DO  0  J  =  2»'3  i./Av^^ftOri 

9  E(J-1»I)=B(J»I)  5  NMGO=C  f  DOT I  K-=l»NOP  WAVE07O0 
IF<NOGO(K)-l ) 12,12,11  WAVE0710 

12  NNGO=NNGO+l                                         '  WAVE0720 
11  CONTINUE  $  IF(NNGO-NOR)  13,14,] 4  ^.VFA710 

13  GO  TO  23  .  '  WAVEO740 

14  CONTINUE  WAVE0.750 
M=l  $  D03.1I  =1  »MM  $  D031J=1,NN  ff  T  F  (  ^CP  (  I  ,  J  )  )  3  1  ,  30  ,  3  1  WAVE0760 

30  XXX  (M)  =  (  1-1  )*GRID  5  YYY  (  v  )  -  (  J-l  )  *OP  I  D  «  M=W|+1  $  \l  =  M-l  WAy€f>77^ 

3  1  CONTINUE  WAVE*7fiO 

CALL  DRAW(N>XXX,YYY,3  .OtLABEL,  T  T  I  T  L  r  ,  10^0.,  10000  .  » 0 »0  ,  2  ♦  2  ,Q,P  ,   V-M/^7^ 

11, LAST)  w-AvFApon 

END  ,.,avcaoio 

SUBROUTINE  RAYN  (  XX  ,  YY  ,  NOP  )  .■  oayN'op?* 
COMMON/ BLK1/X  (  lfiO  )  »Y  (  IOC  )  sCOR^R  (  100  )  .  »HKC  f  lft.0  )  »B  (  3  » 100  )  ,fp<  (  in*  )  ,  PAYNrt5*n 

INOGOt 100) ,NGO,AAA( 100 )  RAYNO840 
COMMON/ 3LK2/T » A»CXY, PAR »3AR »T I VE>GRID».FK ,MAX , A3C ,mv , NN , DEPTH ,KHA ,  PAYMrt85n 

1C0REFA,N0RA,LC,DIST  RAYM0860 

COMMON/ 3LK3/DEP( 100,100) , PDPX , PDPV , PGPX , PCDY , PDDDXX ,PDDPYY  RAYNn/870 

1  ,PDDPXY,PDDPCC ,nDPC  pavmorp.0 

C0MM0lM/SLK4/irFKK  (  100  )  ,  DLO  ,  C  (  6  )  ,'' IT  ,pcL  X,n'cL  v  RAyM^PQO 

IFtf-'AX-l  JXICtllCtlll  PAYNOPon 

110  NGO=2  J  GO  TO  108  daymooio 

111  NGO  =  NOGO  (  NOR  A  )  $  GO  TO  ( 104 ,  1*p  1  NOO  DfivMno;?n 

108  CALL  DEPTFUN(XX»YY)  PAYM0930 
GO  TO  (104,102)NG0  o/\ymaQ4o 

104  NGO  =  1XCCRIF  A  =  0  •  0  53(3   »NCPA ) =0 ,OSFK=C • 0  %    HMA  =  0«C  f  <"-n  Tn  lp3     davmaqcja 
102  IF(MAX-l) 105»105,1C7  PAYM*Q6n 

105  CALL  <FUNCT(A»F<)  5  rO  TO  106  RAYMO970 
107  FK=FF<(NORA)  PAYN^PSO 

106  CALL  '-'OVF  (  XX  »YY  1  1  GO  TOf  1 0&  » !  A9  )  NOO  'ivwrioao 

109  CALL  ?PTA(XX,YV)  o4vM] aaa 
10  3  F FK  (  NOR  A )  -  F  K  I  COR  E  ~  R  (  NOR  A  )  =  G^  °  ^  F  ■-  5  H  H  0  (  N' nD  A-  )  ojvwi^m 

1  ^Hf  |  A  X  NO  ^  0  (  "10  P  A  )  t  "'  ^  ">  '  y  (■''  j«  i  rvv  <r.  y  (  >iosa  i  -vy  <r  AM  i"'n')i  i  r'  0AYM1n2O 

RrTURN  r"*yMi  it> 

j;j r  r?  o 'J  T  I  N  E  ^ "" r  T  F!J  ^!  (  X  X  •  Y  Y  )                      '  ^  c  n  T  ]  n  c  ^ 
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D  I '-'ENS  I  ON  3(6,6)  ,E(6l 


CM  1 A0)  ,D(  100 ) 


C0KMON/BLK1/X  (  100  )  ,  Y  (  1  00  )  ».COR.FcR  (  1  0^  )  ,  "^  (  1<™  )  ,  o  (  ?  , 
1NOGOI ICC) ,NGO»AAA( IOC) 

CCMMON/3LK2/T  »A,CXY, PAR ,3AR »TlME » GRID »FK , MAX  , NOR, MM 
1CCREFA,NORA,'|_0,DIST 

COMMON /3LK 3 /DEP  ( 100 , ICO ) »PDPX ,PD°Y ,PCPX »PCPY »PDDPXX 
1  ,PDDPXY,PDDPCC,PDPC 

COMM0N/BLK4/FFKK(  100)  ,DL0,.C(6)  ,V/,IT  ,DCL  \»DEL  Y 

XR=XX/GRID  I    K=XR+1.5  T  YR=YY/^°!^  «  L=Y°+1»? 

IF( CXX-1.0)*( fMM-1 )-K) )62»7,7 

7  IF(  (YY-1.0)*(  (NN-1  )-l_)  5  62,6,6 
6  D09J=1,6  $  D09I=1,6 

9  D(  I  ,J)=0.0  $  D08  1=1,6 

8  Ed)  =0.0 

M=L-1  S  MA=L+1  $  N=K-1  $  NA=K+1 

DO10  JsMtMA  $  D01C  I=N,NA  $  0 ( I ) = ( I -1 ) *GR I D  $  P(J)= 
D( 1,1)=9.$D(2,1 )=D(2,1  )+0(  I )$D( 3,1 )=D( 3,1 )+P( J)$D(4 
1)*P(J)$D(5,1)=D(5,1)+C(I)**2$D(6,1)=D(6,1)+D(J)**2$ 
2U>J>$Da»2)=D(2»l)$D(2»2)=D(5,l)SD(3»2)=D<4,.l)SD(4 
3)**2*P(J)$D(5»2)=0(5,2)+O(I)**''SD(6»2)=P(6,2)+O(I)* 
4E(2)+0(I)*DEP(I»J)SD(l»3)=P(3,1)$rM2,^)=rM4,l)?rM^, 
53)=D(4,3)+0(I)*P(J)**2SD(t:,3)=rM4,2)SD(6,^)=D(6,^)  + 
6E(3)+P( J)*DEP( I tJ)SO(l ,4)=D(4,l)Sn(2,4)=D(4,2)$0(;?, 
74)=D(4,4)+0(I)**2*P(J)**2SD(5»4)=D(c',4)+O(I)**3*D(j 
8)+0(I)*P(J)**3lE(4)  =  c(4)+0(I)*CMJ)*O^D(I,J)$CMl,^)  = 
9D(5,2)SD(3,5)=D(4,2)SD(4,5)=0(c,4)$D(5,5)=o(5,5)+0( 
1D(4,4)$E(5)=E(5)+0(I )**2*DEP( T,J)$D(1,6)=D(6»1)$D(? 
2,6)=D(6,3)$D(4,6)=D(6,4)$D(5,6)=D(4,4)Sri(6,6)=0(6,6 
3=E(6>+P( J)**2*DEP( I , J ) 

10  CONTINUE 

MPM=7   1  DO  34M=1,6  $  KP  =  0  $  Z  =  0.0  $  D012N=M,6 
1  )  )11,12,12 

11  Z=ABSF(D(N,M) )  S  KP  =  N 

12  CONTINUE  $  IF(M-KP  >13, 20,20 

13  D014J=M,NPM  $Z=D(M,J)  $  D ( M , J ) =D ( KP » J ) 

14  D(KP,J)=Z 
20  IF(AB3F(D(M,M)  ) -.0  0  00 1  )  50  ,  50  ,  3  0 

30  IF(M-6)31,40,40 

31  LP1=M+1  $  D034N  =  LP1,6  $  I F ( D ( N , M ) ) 32 , 34 , 32 

32  RATI0  =  D(N,M)/D(M,M)$  D03  3J  =  l_Pl  ,NPM 

33  D(N,J)=D(N, J)-RATIO*D(M,J) 

34  CONTINUE 

40  00431=1,6  $  11=7-1  $  JPN=7 
143 

41  IIP1=II+1  $  DO  42N=IIP1.6 

42  S=S+D( II ,N)*C(N) 

43  C(  I  I )   =(D(  I  I  ,JPN)-5)/D( I  I ,1  I ) 

50  KER=2 
PRINT  54,XX,YY,KER 

54  FORMAT ( /3X.5HXX  =  1PE20.8 » 3X . 5HYY  =  IPC2'-. 
117H   MATRIX  SINGULAR/) 

51  GO  TO  (52,53)KER 
53  RETURN 


,NN,DE0TH  ,HHA, 
,  Dnnoyv 


( J-l)*ORID 
,  1  )=0(4,  1 )+0( I 
F(  1  )=F( 1  )+D^P 
,2)=P( 4,2 )+0( I 
D(  J)  **?$c(  7.  )  = 
3>=D(6tl) 3^(4, 
P ( J) **3$r( 3  )  = 
4)=0(6,2)?n(4, 
)$D(6»A)=D(6»A 
D(6,1)$D(2,?)= 
I  )**4$D(6»C5  )  = 
»6)=D(6»2  )$D( 3 
)+P( J)**4*F(6) 


$  IMZ-^SFtD^NM) 


S=0.0  $  1^(11-6)41, 43, 


KER=1  5  GO  TO  51 


4X,6H<FD  =  12 


riCDTl'>7A 

DFPTln80 

DEDTll^o 
DpPTl 1 10 
Dcn>T1120 
Dpdti  1  30 
D^DJl 1 4^ 
PpPTl  1?0 
DFPT1160 
DFPT1170 
DFPT1 1RO 
DFPT11PO 
DEDT1200 
DEPT1 210 
DEPT1220 
DFDTI 230 
0PPT124O 
DPPTl ?K0 
op°T1 26° 
Dpdti 27o 
PicDj]  ?PO 
DCDT1 ?P0 

0FDJ1  -JOO 

nPDJl^lO 
DEPT1 320 
nc-p-n  »30 

DEPT1340 
DroT1^50 
OPPT1360 
DPPT1370 
HFPT1 330 

DFPT1 4*>0 

0PPT1410 
DF°T1420 
DFPT1430 
DFDT1440 
DCPT1450 
DFPT1460 
DEPT1470 
DcpT14«0 
DEPT1490 
DFPT1500 
PFPTISIO 
DFPT1 520 
0FPT153O 
DFOT1R40 
DEPT15*o 
DEDT1560 
DEPT1570 
DFPT158'> 
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62 

NG0  =  1  $  GO  TO  81 

64 

CONTINUE 

81 

RETURN 

END 

SUBROUTINE  KFUNCT 

52  CONTINUE  DFPT1590 

DEPTH=C(  1  )+C(2  )  *XX+C  (  3  )*YY+C(4]  *XX*YY+C(  c  )*XX**2+C(  6)*YY**2        D^PTl  600 
IF(DEPTH)62,62,65  DrPT1610, 

65  AL0=PAR*T  X  DLO=DEPTH/ ALO  $  I  F  (.DLO-0  .  5  )  63  »63  *66  DFPT1620 

66  CXY=PAR  J  L0=1  I  NGO  =  2  I    FK=0.  $  GO  TO  67  DEPT1630 
63  DO60M=1.5C  I  CXX=PAR*TANHF ( BA^*DEPTK/CXY )                           DFPT164H 

IF( ABSF(CXX-CXY)-0.01 ) 61 .61  ,60  DFPT1650 

60  CXY=(CXX+CXY)*.5  D-cPT]66n 

61  CONTINUE   $  NG0=2  1  L0=2  DEPT1670 

67  PDPX=C( 2)+C (4)*YY+2*C( 5)*XX  DEPT1680 
PDPY=C(3)+C(4)*XX+2*C(6)*YY  -  DEPT1690 
PDDPXY=C(4)  $  PDDPXX=2*C( 5 )  $  PDDPYY=2*C ( 6 )  DFPT1700 
GO  TO  64            ■   v  DC.^T1710 

DFPT1720 

DEPT1730 

DFojl74n 

DEPT1750 

( A,FK  )  KFUN176* 

COMMON/ 3LK1/X(  IDO')  »Y  (  IOC)  »CCRFFR(  100  )  »HHO(  100  )  » B  (  3  ,  100  )  ,FFK(  100  )  i  KFUN1770' 

lNOGO(lOO)  ,NGO,AAA(  IOC ■  J  K FUN  1780 

C0MM0N/3LK2/T,UU»CXY,PAR,3ARrT  TIE, GRID, VV .MAX .NOR .MM » NN . DEPTH ,HHA  K FUN 1790 

1  ,COREFA»NORA,LO,DIST  KFUN1800 

COMMON/ BLK3/DEP( 100. 100).  .PDPX , PDPY , PCPX , PCPY , PDDPXX , PDDPYY         K FUN 1810 

l.PDDPXY.PDDPCCPDPC  KFUN1820 

3  GO  TO  ( 5,6)L0  _  KFUM1830 
6  Rl  =  CXY/32.2  $  R2  =  R1**3*BAR'**2  •  $  R3=Rl**5*BAR**4  J  94  =R  l**7.*PAP**6  KFUN1840 

PDPC=2.*R1+4.*R2/3.+6.*R3/5.+p.*R4/7.  '      KPUN1850 

PDDPCC  =  <2.*R1M,*R2+6.*R3  +  8.*R4) /CXY  K FUN I8  60 

PCPX=PDPX/PDPC  J  PCPY=PDPY/PDPr  KFUN1870 

FK  = (PCPX*SINF( A)-PCPY*COSF(A ) ) /CXY  I  GO  TO  4  KFUN1880 

5  FK=0.  KFUN1890 

4  RETURN  KFUN1900 
END  KFUN1910 
SUBROUTINE  MOVE(X.Y)  MOVE1920 
C0MM0N/BLK1/U( 100 ) ,V( 100) ,CORrFR( 100  )  ,HHO  (.100  )  »B  (  3.100  )  »FFK(  100  )  »  MOVE  19^0 

1N0G0( 100) ,NGO.AAA( 100)  MOVE1940 

C0MM0N/BLK2/T.A.CXY.PAR.BAR.TIME.GRID.FK.MAX , NOR , MM , NN , DEPTH ,HHA ,  MOVE  19  50 

1COREFA,NORA,LO.DIST  M0VF196n 

COMMON/ BLK3/DEP( IOC, 10  0) , PDPX , PDPY , PCPX , PCPY , ^DDPXX » PPDPYY         ^0VC1 07^ 

l.PDDPXY.PDDPCCPDPC  ^OVF1<?80 

C0MM0N/BLK4/FFKK( 100) ,DLO, C(6) , "IT, PEL  X,DFL  Y  ^OVF1990 

FKBAR  =  FFKK(NCRA)  MOVF'2^^0 

IF(MAX-l) 1.1.4  MOVE2010 

1  FKBAR=FK  MnvF2^20 

4  MIT=1  MOVE2030 

DEL  D=TIME*CXY  MOVE2040 

GO  TO  (22.2DL0  MOVE2O50 

22  XX=X+DELD*COSF( A)  5  YY=Y+DELP*S I NF ( A )  $  AA=A  $  FK<=0.  $  FKBAP  =  C   wOVr2n60 

GO  TO  6  M0VE2O7O 

21  DO  20IT=1 ,15  ^OVc2^80 

19  DEL  A  =  FKBAR*DEL  D  1    £A=A+DEL  ft  <E  ABAR  =  A  + .  5*DEL  A  S  DFL  X=PCL  D*    MOVc209O 

1  COSF(ABAR)  $  DrL  Y  =  ~FL  0*S  I N1-  (  A  BA  9  )  $  X.X  =  X+,OFL  X  $  YY  =  v  +  ^cL  Y     M0vr2]*A 

GO  TO  (10 1,6)  MIT  •  ^OVE2110 
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% 


101  CALL  DEPTFUN(XX.YY)      $  GO  TO ( 38 , 1 0 ) NGO  MOVE2120 

10  CALL  KFUNCT(AA,FKK )  $  FKEAR= . 5* < FK+FKK >  MOVE2130 

IF( IT-13)5i37.9  MOVF2140 

37  FKKPP=FKBAR  MOVE2150 

5  IF(MAX-1)7,7,9  MOVE216D 

7  IF( IT-1)20,20,9  MOVE2170 
9  IF(ABSFtFKKP-FKBAR)-. 00001)6, 6, 20  MOVF2180 

20  FKKP=FK8AR  MOVE2190 

24  IF(ABSF(FKKPP-FK8AR)-. 00001)18, 18, 17  MOVE2200 

17  MIT=3  $  NGO=l  $  GO  T038  MOVF2210 

18  FKBAR=.5*(FKBAR+FKKP)  £  MIT=2  $  GO  TO  19  MOVF2220 

6  NGO=2  $  GO  TO  8  MOVE7230 

8  X=XX  $  Y  =  YY  $  A  =  AA  $  FK=FKK  M0VE224* 

38  CONTINUE  M0VE2250 
FFKK(NORA)=FKBAR  ^  MOVE2260 
RETURN  MOVE2270 
END  MOVE2280 

'SUBROUTINE  BETA(XX.YY)  BETA2290 

COMMON/BLK1/X(100) »Y( 100) »CORPFR( 100),HH0( 100 ) »B ( 3 » 100 ) ,FFK ( 100),  BET A2 300 

.  INOGOdOO)  »NGO»AAA(  100)  BET A2 3 10 

-:?  COMMON/  BLK2/T»A,CXY,PAR,BAR,TTMF»GRID»FK,MAX,N0R,  VIM,  NN,  DEPTH, HM  A,  BET A2 3 20 

1COREFA,NORA,LO,DIST                   _  BETA2330 

COMMON/ BLK3/DEP(  100,100)  ,PDPX  ,  PDPY  ,PCPX  ,PCPY  ,PDDPXX.,PDDPYY  BETA2  340 

lfPDDPXY.PDDPTCPDPC  BETA2350 

C0MM0N/BLK4/FFKK( 100) ,DL0»C(6) ,MlT,DEL  X,DEL  Y  BETA2360 

GO  TO  (5»6)L0  BETA2370 

5  COREFA  =  l.  $  HHA  =  1.  $  B  (  3  ,N0RA  )  =B  (  2  ,NORA  )  I  COR"EFB  =  l.  $  HHB=1.  BEJA2380 
GO  TO  7  BETA.2390 

6  PCCPXX=PDDPXX/(PDPC+PDDPCC)  BETA2400 
PCCPYY*PDDPYY/(PDPC+PDDPCC)  BETA2410 
PCCPXY»PDDPXY/(PDPC+PDDPCC)  BETA2420 
P=(-COSF( A)*PCPX-SINF(A)*PCPY) /CXY   $  0= ( S I NF ( A ) **2*PCCPXX-2 .*     BETA2430 

15INF(A)*C0SF(A)*PCCPXY+C05F(A)**2*PCCPYY)/CXY  BET A? 440 

DD=SQRTF(  (DEL  X)**2+(DEL  Y)**?)  BETA2450 

Bl3»N0RA)=(B(l»N0RA)*(P*DD-2. ) +B ( 2 »NORA ) * ( 4.-2.*DD**2*Q ) )  BETA2460 

l/(2*+P«DD)  BETA2470 

C0REFA»1./SQRTF(ABSF(B(2»N0RA) ) )  BETA2480 
CC0*CXY/PAR  $  HSH0L=3. 2519-12. 8l50*CC0+28. 81 12*CC0**2-29. 9257*CC0  BETA2490 

1*»3+11.6815*CC0**4  $  HHA=C0REFA  *  HSHOL  BETA2500 

7  RETURN  BETA2510 
END  BETA2520 
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-9 
-8 

7 
-6 


Example  of  Input 

X         Y     Bl   B2   T   A    NOR    OIST    TIME   GRID      MM    NN 
4500.     33000.   1.  1.  20.-35.   18     1500.    30.   1500.      52   56 

440.  460.  515.  505.  510.  470.  425.  355.  300.  330.  370.  360.  340.  330.  -11 

340.  345.  350.  365.  375.  390.  455.  560.  700.  850.  950.  865.  705.  630.  -10 

550.  490.  470.  460.  415.  380.  335.  400.  460.  520.  560.  595.  620.  590. 

510.  460.  400.  315.  275.  300.  220.  140.   80.   55.   53.   51.   51.   51. 

360.  405.  450.  435.  445.  380.  330.  270.  235.  250.  280.  280.  260.  270.  -7 

290.  297.  300.  305.  315.  360.  420.  490.  655.  810.  930.  870.  750.  620.  -o 

540.  500.  515.  515.  480.  445.  440.  500.  540.  593.  630.  630.  600.  590.  -5 

610.  560.  515.  400.  315.  275.  234.  190.  110.   63.   53.   50.   50.   50.  -4 

295.  370.  370.  325.  325.  307.  260.  190.  170.  175.  190.  190.  180.  205.  -3 

240.  250.  250.  265.  275.  305.  405.  470.  570.  750.  900.  885.  800.  630.  -2 

560.  540.  550.  560.  550.  540.  550.  605.  620.  640.  610.  570.  530.  510.  -1 

540.  590.  570.  430.  300.  310.  280.  240.  170.   70.   55.  50.   50.  50.  0 

240.  250.  235.  195.  225.  200.  190.  160.  130.  140.  160.  170.  170.  180.  1 

190.  195.  210.  230.  250.  280.  350.  427.  490.  670.  820.  860.  820.  680.  2 

610.  590.  590.  615.  600.  610.  620.  660.  650.  620.  560.  505.  480.  470.  3 

480.  520.  560.  515.  405.  370.  320.  260.  190.   80.   60.   49.   49.   49.  4 

150.  195.  200.  110.  125.  120.  115.   95.   80.   98.  105.  110.  115.  135.  5 

162.  170.  180.  205.  210.  210.  290.  380.  455.  580.  740.  785.  840.  720.  6 

660.  650.  640.  660.  660.  670.  690.  685.  650.  585.  510.  480.  465.  440.  7 

450.  490.  550.  525.  450.  380.  325.  270.  215.  120.   65.   50.   50.   50.  8 

80.   80.   90.   75.   73.   70.   72.   70.   70.   70.   72.   80.   85.   95.  9 

110.  120.  140.  160.  175.  195.  240.  290.  390.  475.  610.  670.  750.  815.  10 

760.  720.  715.  725.  730.  710.  680.  655.  620.  540.  480.  445.  430.  420.  11 

440.  410.  550.  490.  415.  370.  325.  285.  235.  170.  85.   65.   50.   50.  12 

61.   60.   60.   58.   60.   60.   61.   63.   60.   60.   65.   65.   75.   75.  13 

80.   85.   100.  120.  130.  120.  160.  230.  280.  380.  460.  540.  630.  775,  14 

770.  750.  735.  715.  665.  635.  620.  610.  570.  485.  440.  415.  385.  375.  15 

425.  490.  550.  500.  445.  375.  315.  295.  245.  190.  110.  80.   65.   50.  16 

47.   47.   48.   49.   52.   54.   53.   54.   56.   59.   60.   60.   60.   62.  17 

63.   65.   70.   77.   82.   85.   110.  150.  185.  250.  340.  395.  465.  550.  18 

560.  610.  590.  600.  570.  540.  525.  520.  500.  460.  410.  380.  360.  365.  19 

420.  510.  515.  500.  465.  410.  283.  265.  235.  210.  145.  105.  60.   52.  20 

38.   38.   38.   41.   46.   48.   51.   53.   54.   56.   57.   57.   59.   60.  21 

60.   62.   63.   67.   70.   70.   80.   95.  115.  160.  210.  300.  370.  450.  22 

450.  465.  490.  510.  470.  465.  460.  460.  440.  415.  395.  385.  370.  360.  23 

430.  520.  470.  430.  415.  420.  240.  270.  254.  260.  250.  215.  160.   55.  24 

31.   32.   33.  38.    40.   43.   48.   50.   51.   54.   55.   56.   57.   58.  25 

59.   60.   61.  62.    63.   64.   65.   70.   80.   86.  140.  185.  285.  385.  26 

360.  380.  420.415.   400.  390.  390.  410.  380.  365.  350.  350.  350.  352.  27 

500.  470.  410.360.   350.  370.  380.  375.  340.  310.  315.  300.  195.  120.  28 

30.   28.   33.  36.    35.   35.   38.   43.   46.   49.   51.   53.   55.   56.  29 

58.   59.   60.  60.    62.   62.   63.   64.   65.   73.   80.  120.  185.  260.  30 

240.  280.  325.  325.  325.  330.  340.  330.  310.  290.  305.  320.  405.  450.  31 

460.  410.  300.  260.  260.  295.  300.  290.  300.  300.  290.  310.  300.  220.  32 

13.   24.   26.   28.   28.   31.   34.   38.   42.   45.   48.   51.   54.   54.  33 

56.   57.   58.   59.   59.   60.   60.   62.   64.   66.   70.  120.  160.  160.  34 

155.  180.  240.  240.  230.  230.  260.  270.  270.  280.  325.  370.  440.  470.  35 

420.  370.  300.  240.  240.  210.  180.  200.  200.  190.  210.  220.  300.  260.  36 

18.   20.   24.   25.   25.   26.   27.   33.   36.   40.   44.   49.   51.   53.  37 

55.   56.   58.   58.   57.   59.   60.   60.   62.   63.   63.   85.  120.  105.  38 
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Example  of   Output 


XV 

YV           Bl    82 

T 

Al 

NOR 

OIST           TIME 

GRID         MM      NN 

U 500.0 

33000.0      1       1 

20 

-35 

18 

1500.00    30.0 

1500.0      52      56 

NUMBER    Of    CREST    ■ 


X 

Y 

COREFR 

HHO 

NGO 

X 

Y 

2339.60 

10412.11 

1.000 

1.000 

2 

2626.38 

10821.68 

2913.17 

11231.26 

1.000 

1.000 

2 

3199.96 

11640.84 

3486.75 

12050.41 

1.000 

1.000 

2 

3773.54 

12459. S9 

4060.33 

12869.56 

1.000 

1.000 

2 

U347.ll 

13279.14 

4633.90 

13688.72 

1.000 

1.000 

2 

U920.69 

14098.29 

5207.48 

14507.87 

1  .000 

1.000 

2 

5494.27 

14917.44 

5781.05 

15327.02 

1.000 

1.000 

2 

6067.84 

15736.60 

6354.63 

16146.17 

1.000 

1.000 

2 

6641.42 

16555.75 

6928.21 

16965.33 

1.000 

1.000 

2 

7215.00 

17374.90 

COREFR  HHO 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 


NGO 
2 
2 
2 
2 
2 
2 
2 
2 
2 


NUMBER    OF    CREST    - 


X 

Y 

COREFR 

HHO 

NGO 

X 

Y 

3179.19 

982*. 22 

1.000 

1.000 

2 

3465.98 

10233.79 

3752.77 

10643.37 

1.000 

1.000 

2 

4039.56 

1 1052.94 

4326.35 

11462.52 

1.000 

1.000 

2 

4613.13 

11872.10 

4899.92 

12281.67 

1.000 

1.000 

2 

5186.71 

12691.25 

5473. SO 

13100.83 

1.000 

1.000 

2 

5760.29 

13510.40 

6047.07 

13919.98 

1.000 

1.000 

2 

6333.86 

14329.55 

6620.65 

14739.13 

1.000 

1.000 

2 

6913.46 

15164.19 

7182.28 

15566.65 

1.000 

.899 

2 

7454.38 

15986.51 

7728.66 

16404.84 

1.000 

.895 

2 

8001.84 

16823.95 

COREFR  HHO 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  1.000 

1.000  .905 

1.000  .897 

1.000  .895 


NGO 
2 
2 
2 
2 
2 
2 
2 
2 
2 


NUMBER  OF   CREST 


X 

Y 

COREFR 

HHO 

NGO 

X 

Y 

COREFR 

HHO 

NGO 

4005.32 

9245.76 

1.000 

.902 

2 

4292.80 

9654.85 

1.000 

.899 

2 

4577.T1 

10065.74 

1.000 

.897 

2 

4858.41 

10479.58 

1.000 

.897 

2 

5140.28 

10890.72 

1.000 

.897 

2 

5423.42 

11304.73 

1.000 

.898 

2 

5711.14 

11713.66 

1.000 

.898 

2 

6001.92 

12096.33 

1.000 

.898 

2 

6307.35 

12513.13 

1.000 

.903 

2 

6587.17 

12931.41 

1.000 

.903 

2 

6854.59 

13354.55 

1.000 

.898 

2 

7124.80 

13775.74 

1.000 

.896 

2 

7394.65 

14197.17 

1.000 

.895 

2 

7694.29 

14666.16 

1.024 

.917 

2 

7930.90 

15056.20 

1.004 

.899 

2 

8191.41 

15480.48 

1.032 

.923 

2 

8446.16 

15904.45 

1.031 

.922 

2 

8727.76 

16319.95 

1.009 

.903 

2 
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NUMBER  OF  CREST  -   22 


X 

Y 

COREFR 

hhO 

NGO 

X 

Y 

COREFR 

HHO 

NGC 

10264.28 

4616.24 

0 

0 

1 

10853.02 

4550.71 

0 

0 

1 

11438.10 

4352.73 

0 

0 

1 

12072.87 

3316.63 

0 

0 

1 

12920.26 

2737.33 

0 

0 

1 

15666.05 

1290.81 

.270 

.291 

2 

17781.43 

2925.94 

.467 

.486 

2 

16546.11 

1779.73 

8.942 

9.377 

2 

18531.25 

4023.30 

.882 

.924 

2 

18968.07 

4701.04 

.769 

.805 

2 

19101.90 

4934.51 

.757 

.791 

2 

19511.47 

5662. 92 

.859 

.900 

2 

19789.68 

6181.58 

1.139 

1.199 

2 

20449.56 

7333.27 

1.204 

1.271 

2 

20588.90 

7626.71 

.930 

.986 

2 

20819.12 

8054.09 

2.473 

2.659 

2 

20812.04 

8068.66 

1.437 

1.542 

2 

21047.31 

8517.89 

1.045 

1.125 

2 

NUMBER  OF  CREST  •   23 


X 

Y 

COREFR 

HHO 

NGO 

X 

Y 

COREFR 

HHO 

NGO 

10264.28 

4616.24 

0 

0 

1 

10853.02 

4550.71 

0 

0 

1 

11438.10 

4352.73 

0 

0 

1 

12072.87 

3316.63 

0 

0 

1 

12920.26 

2737.33 

0 

0 

1 

15818.45 

860.81 

.252 

.287 

2 

18177.11 

2611.34 

.438 

.488 

2 

16808.11 

1362.17 

6.696 

7.344 

2 

18949.52 

3760.33' 

.682 

.962 

2 

19386.59 

4438.11 

.756 

.840 

2 

19523.12 

4669.33 

.725 

.807 

2 

19940.13 

5417.85 

.850 

.949 

2 

20218.44 

5951.79 

1.143 

1.270 

2 

20883.05 

7120.80 

1.131 

1.285 

2 

21014.87 

7416.37 

.942 

1.068 

2 

21227.05 

7848.73 

2.254 

2.625 

2 

21212.27 

7844.86 

1.365 

1.584 

2 

21454.05 

8311.97 

1.031 

1.228 

2 

NUMBER  OF  CREST 


24 


X 

Y 

COREFR 

HHO 

NGO 

X 

Y 

COREFR 

HHO 

NGO 

10264.28 

4616.24 

0 

0 

1 

10853.02 

4550.71 

0 

0 

1 

11438.10 

4352.73 

0 

0 

1 

12072.87 

3316.63 

0 

0 

1 

12920.26 

2737.33 

0 

0 

1 

15936.15 

482.73 

.238 

.296 

2 

18502.09 

2350.98 

.415 

.536 

2 

17039.75 

997.51 

5.649 

7.103 

2 

19320.29 

3526.00 

.863 

1.084 

2 

19743.52 

4220.11 

.753 

.964 

2 

19879.31 

4454.54 

.708 

.915 

2 

20297.47 

5212.49 

.843 

1.131 

2 

20591.71 

5763.85 

1.124 

1.408 

2 

21240.91 

6950.72 

1.088 

1.493 

2 

21375.95 

7248.08 

.945 

1.397 

2 

21569.75 

7698.52 

2.148 

3.043 

2 

21548.37 

7673.71 

1.332 

1.864 

2 

21784.38 

8176.83 

1.041 

1.774 

2 

36 


APPENDIX  II 


SUBROUTINE  DRAW 


C 

C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


(NUMPTS,  X,  Y,  MODCURV,  ITYPE,  LABEL, 
ITITLE*  EXSCALE*  YSCALE.  IXUP,  IYRIGHT, 
MODEXAX.  MODEYAX,  IWIDE.  IHIGH»  IGRID, 
LAST) 


A  GENERAL  CURVE 
PROGRAMMER 
DATE 
SYSTEM 
OUTPUT 
NOTE 

INPUT  ARGUMENTS 

1,  NUMPTS 


2.  X 


DRAWING  AND  POINT  PLOTTING  SUBROUTINE 
J.  R.  WARD 
FEB.  1964, 
FORTRAN  60 

LOGICAL  TAPE  NUMBER  8 
ASTERISKS  MARK  CHANGES 


REVISED  JUNE  1965 


FOR  FORTRAN  63 


3.  Y 


4.  MODCURV 


5.  ITYPE 


6.  LABEL 


AT 
IN 


LEAST  EQUAL 
CALLING 


OOOOOOOO 

0OOOOO10 

00000020 

30 

0     40 

50 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

000   120 

00000130 

000   140 

NUMBER  OF  POINTS  TO  BE  PLOTTED,  THIS  MUST  ALWAYS00000150 

BE  AT  LEAST  2»  AND  MUST  NOT  EXCEED  30  FOR  POINT  00000160 

PLOTTING,  OR  900  FOR  CURVE  DRAWING,  00000170 

000  180 
00000190 
00000200 
00000  210 
000  220 
00000230 
00000240 
000  250 
00000260 
00000270 
00000280 
00000290 
GRAPH00000300 
00000310 
000  320 
00000330 
POINTS00000  340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000400 
,00000410 
00000420 


ARRAY  OF  X-ORDINATES,  DIMENSION 
TO  NUMPTS  AND  NOT  MORE  THAN  900 
PROGRAM. 


ARRAY  OF  Y-ORDINATES.  DIMENSION  AS  FOR  THE 
X  ARRAY  IN  THE  CALLING  PROGRAM. 

CONTROLS  THE  NUMBER  OF  CURVES,  AND/OR  SETS 
POINTS,  ON  EACH  GRAPH.  THE  CODES  ARE 


OF 


0 
1 
2 
3 

CONTROLS 
0 


ONLY  ONE  PLOT  ON  THIS  GRAPH 
FIRST  PLOT  ON  MULT  I -PLOT  GRAPH 
INTERMEDIATE  PLOT  ON  MULTI-PLOT 
LAST  PLOT  ON  MULT  I -PLOT  GRAPH. 

THE  TYPE  OF  PLOT.  THE  CODES  ARE 
STRAIGHT  LINES  JOIN  SUCCESSIVE 
(STANDARD  CURVE  DRAWING) 
POINTS  PLOTTED  WITH  CROSS  (X) 

PLOTTED 

PLOTTED 

PLOTTED 

PLOTTED 


POINTS 
POINTS 
POINTS 
POINTS 


WITH 
WITH 
WITH 
WITH 


PLUS  (+) 
SQUARE 
DIAMOND 
TRIANGLE 
( ITYPE»1  THRU. 


5) 


1 

2 

3 

4 

5 
WHEN  POINTS  ARE  BEING  PLOTTED 
THE  POINTS  ARE  NOT  CONNECTED. 

000  430 
IF  A  CURVE  IS  BEING  DRAWN  (ITYPE  ■  0),  LABEL  IS  00000440 
A  4-CHARACTER  BCD  CURVE  IDENTIFIER  WHICH  WILL  BE00000450 
REPRODUCED  AT  THE  END  OF  THE  CURVE.  FOR  EX AMPLE ,00000 460 
LABEL  =  4H  ONE  •  IF  POINTS  ARE  BEING  PLOTTED,  00000470 
LABEL  IS  AN  8-CHARACTER  IDENTIFIER.  THE  FIRST  4  00000480 
CHARACTERS  ARE  REPRODUCED  WITH  THE  FIRST  PLOT TEDOOOOO 490 
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c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


POINT,  AND  THE  LAST  4  WITH  THE  LAST  POINT,  SET 
TO  BLANK  ANY  UNWANTED  CHARACTERS. 


7.  ITITLE 


8.  EXSCALE 

9.  YSCALE 

10.  IXUP 

11.  IYRIGHT 

12.  MODEXAX 

13.  MODEYAX 


14.  IWIDE 

15.  IHIGH 

16.  IGRID 

17.  LAST 


00000500 
00000510 
000  520 
00000530 
00000540 
00000550 


AN  ARRAY  OF  TWELVE  8-CHARACTER  BCD  WORDS. 
THE  FIRST  SIX  WORDS  WILL  BE  REPRODUCED  AS  THE 
FIRST  LINE  OF  GRAPH  TITLE.  AND  THE  LAST  SIX 
WORDS  WILL  FORM  THE  SECOND  LINE.  THE  TITLE  MUST  00000560 
INCLUDE  THE  USERS  JOB  IDENTIFICATION.  DIMENSION  00000570 
12  IN  CALLING  PROGRAM.  AND  SET  TO  BLANK  ALL  00000580 
UNWANTED  CHARACTERS.  00000590 

000  600 
X-SCALE  (UNITS/ INCH)  AS  POSITIVE  FLOATING  POINT  00000610 
VARIABLE  WITH  ONE  FIGURE  SIGNIFICANCE.  SET  TO  00000620 
ZERO  FOR  AUTO-SCALE.  00000630 

000  640 
Y-SCALE  (UNITS/INCH)  AS  POSITIVE  FLOATING  POINT  00000650 
VARIABLE  WITH  ONE  FIGURE  SIGNIFICANCE.  SET  TO  00000660 
ZERO  FOR  AUTO-SCALE.  00000670 

000  680 
00000690 
00000700 
00000710 
000  720 
Y-AXIS  OFFSET  FROM  LEFT  OF  GRAPH  IN  INCHES.  THIS00000730 
MUST  NOT  EXCEED  IWIDE »  AND  MUST  NOT  BE  NEGATIVE. 00000 740 

000  750 
MODE  OF  X-AXIS  OFFSET.  SEE  CODES  BELOW.  00000760 

000  770 
MODE  OF  Y-AXIS  OFFSET.  THE  CODES  ARE  AS  FOLLOWS  00000780 

0  COMPUTED  OFFSET.  HOLDING  ORIGIN  ON  00000790 
GRAPH.  THE  CORRESPONDING  IXUP  OR  00000800 
IYRIGHT  IS  IGNORED  00000810 

1  COMPUTED  OFFSET.  WITH  ORIGIN  OFF  THE  00000820 
GRAPH  IF  APPROPRIATE.  THE  CORRESPOND-00000830 
ING  IXUP  OR  IYRIGHT  IS  IGNORED.  USE  00000840 
ONLY  WITH  AUTO-SCALE 

2  AXIS  OFFSET  AS  SPECIFIED  BY  IXUP  OR 
IYRIGHT. 


X-AXIS  OFFSET  FROM  BOTTOM  OF  GRAPH  IN  INCHES. 
THIS  MUST  NOT  EXCEED  IHIGH.  AND  MUST  NOT  BE 
NEGATIVE. 


00000850 
00000860 
00000870 
000  880 
00000890 
00000900 
000  910 
HEIGHT  OF  GRAPH  IN  INCHES.  THIS  MUST  NOT  EXCEED  00000920 


WIDTH  OF  GRAPH  IN  INCHES.  THIS  MUST  NOT  EXCEED 
NINE.   ZERO  WILL  BE  READ  AS  EIGHT  INCHES. 


FIFTEEN.  ZERO  WILL  BE  READ  AS  EIGHT  INCHES. 

IF  SET  TO  1,  A  ONE  INCH  BY  ONE  INCH  GRID  WILL 
BE  SUPERIMPOSED  ON  THE  GRAPH. 

INDICATES  TO  CALLING  PROGRAM  WHETHER  LAST  PLOT 
WAS  COMPLETED  SUCCESSFULLY.  THE  CODES  ARE 

0  LAST  PLOT  COMPLETED  SUCCESSFULLY 

1  LAST  PLOT  NOT  SUCCESSFUL 

2  LAST  PLOT  NOT  SUCCESSFUL  AND  NO 


00000930 
000  940 
00000950 
00000960 
000  970 
00000980 
00000990 
00001000 
00001010 
00001020 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


FURTHER  GRAPH  OUTPUT  WILL  BE  ATTEMPT- 
ED UNTIL  MODCURV  IS  NEXT  ONE  OR  ZERO 
3     DRAW  WAS  ENTERED  WITH  MODCURV  NOT 

EQUAL  TO  ONE  OR  ZERO  WHILE  THE  ER«OR 
LOCK-OUT  WAS  SET. 

THIS  ARGUMENT  MUST  ALWAYS  BE  A  NAME  IN  THE  CALL 

STATEMENT.  NEVER  A  NUMBER. 


NOTE 

ALL  ARGUMENTS  FROM  NUMBER  7  THRU.  NUMBER  16  ARE  IGNORED  WHEN 
MODCURV  IS  EITHER  2  OR  3.  HOWEVER*  ARGUMENTS  MUST  NEVER  BE 
OMITTED  FROM  THE  CALLING  STATEMENT.  IT  IS  MERELY  THEIR  VALUES 
WHICH  ARE  THEN  IRRELEVANT.  ARGUMENTS  MAY  BE  LISTED  BY  NAME  OR 
VALUE  IN  THE  CALL  STATEMENT.  NO  VALUE  IN  THE  CALLING  PROGRAM 
WILL  BE  ALTERED  BY  THIS  SUBROUTINE. 

REFERENCE  

THE  BINARY  TAPE  FORMAT  REQUIRED  BY  THE  OFF-LINE  PLOTTER  IS 
DESCRIBED  IN  THE  WRITEUP  OF  THE  CDC  160A  GRAPH  PLOT  PROGRAM 
(IDENT.  B001).  THE  FORMAT  REQUIRED  BY  THE  CDC  160  PROGRAM  IS 
SIMILAR  EXCEPT  THAT  THE  INTERPOLATION  ARGUMENT  MUST  BE  ZERO. 


* 


♦ 
C 
C 

c 
c 
c 


DIMENSION  X(900),  Y(900).  lTITLEC12)t  JXTIT(12)t  JYTITU2). 

1  LTITLEU4),  KAXISC5),  ICURVC460),  JGRID(25)»  ICONT(l), 

2  JJTITLEU2) 
IPOINT  =  ITYPE 

CON(  ICONTRL  =  40000B). 

CON(ICURV3  =  3777377720202020B.  ICURV4  »  0104000000000000B) • 

REPLACE  WITH  DATA  STATEMENT  IN  FORTRAN  62-3. 
PUT  ITEST  =  0  IN  DATA  STATEMENT. 
IF  (JTEST  -  73546912)   9070 ,9071 ,9070 

9070  ITEST  =  0 
JTEST  »  73546912 

9071  CONTINUE 
REMOVE  ABOVE  NONSENSE  IN  FORTRAN  63. 

CHECK  PREVIOUS  OPERATION  OF  ROUTINE,  IF  ANY.  CODES  ARE 
ITEST  =0     IF  PREVIOUS  GRAPH,  IF  ANY,  COMPLETED 
ITEST  =1     IF  PREVIOUS  GRAPH  NOT  COMPLETED 
ITEST  =2     IF  ERROR  FOUND  WHILE  MODCURV  WAS  ONE,  OR  IF 
MODCURV  WAS  ILLEGAL. 
IFUTEST  -  2)1000,1001,1000 

1001  IF(MODCURV) 1003,1002,1003 

1002  ITEST  x  o 
GO  TO  1000 

1003  IF(MODCURV  -  1)1004,1002,1004 

1004  LAST  ■  3 
RETURN 

SET  UP  ERROR  RETURN  ROUTINE.  ENTRY  AT  STATEMENT  1005. 

1005  IF( ITESTJ1009, 1006, 1009 

1006  IF(MODCURV) 1007,1008,1007 

1007  PRINT  1100 


00001030 
00001040 
00001050 
00001060 
00001070 
00001080 
00001090 
00001100 
00001110 
00001120 
00001130 
00001140 
00001150 
00001160 
00001170 
00001180 
00001190 
00001200 
00001210 
00001220 
00001230 
00001240 
00001250 
00001260 
00001270 
00001280 
00001290 
00001300 
00001310 
00001320 
00001330 
00001340 
00001350 
00001360 
00001370 
00001380 
00001390 
00001400 
00001410 
00001420 
00001430 
00001440 
00001450 
00001460 
00001470 
00001480 
00001490 
00001500 
00001510 
00001520 
00001530 
00001540 
00001550 
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1100 

FORMAT  (  59H  NO  FURTHER  GRAPH  OUTPUT  UNTIL  MODCURV  NEXT  IS  ZERO  OR0O001560 

] 

L  ONE.  »/) 

00001570 

ITEST  =  2 

00001580 

LAST  «  2 

00001590 

RETURN 

00001600 

1008 

PRINT  1101 

00001610 

1101 

FORMAT  (  30H  THIS  PLOT  WILL  NOT  BE  OUTPUT,  ./)                      00001620 

LAST  »  1 

00001630 

RETURN 

00001640 

1009 

IF(MODCURV  -  2) 1010, 1008*1010                                       00001650 

1010 

IFtMODCURV  -  3)1007,1011,1007                                       00001660 

1011 

ITEST  =  0 

00001670 

GO  TO  1008 

00001680 

CHECK  LEGALITY  OF  INPUT  ARGUMENTS.                              00001690 

1000 

IF(NUMPTS  -  2)1.2.2 

00001700 

1 

PRINT  100 

00001710 

100 

FORMAT  (/,  32H  NUMPTS  MUST 

NOT  BE  LESS  THAN  2.  )                    00001720 

GO  TO  1005 

00001730 

2 

IF( IPOINT)9000. 9004. 9001 

00001740 

9000 

PRINT  9100 

00001750 

9100 

FORMAT  </»  15H  ILLEGAL  ITYPE.  )                                      00001760 

GO  TO  1005 

00001770 

9001 

IFUPOINT  -  5)9002.9002.9000                                         00001780 

9002 

IFCNUMPTS  -  30)3.3,9003 

00001790 

9003 

PRINT  9101 

00001800 

9101 

FORMAT  (/.  46H  NUMPTS  MUST 

NOT  EXCEED  30  FOR  POINT  PLOTTING.  )     00001810 

GO  TO  1005 

00001820 

9004 

IF(NUMPTS  -  900)3.3.9005 

00001830 

9005 

PRINT  9102 

00001840 

910  2 

FORMAT  (/,  28H  NUMPTS  MUST 

NOT  EXCEED  900.  )                        00001850 

GO  TO  1005 

00001860 

3 

IX  =  1HX 

00001870 

IY  =  1HY 

00001880 

AMAXX  »  -0.2E+100 

0000 1 890 

AMAXY  «  -0.2E+100 

00001900 

AMINX  »  +0.2E+100 

oooomo 

AMINY  ■  +0.2E+100 

00001920 

DO  1020  I«  1 .NUMPTS 

00001930 

AMAXX  »  MAX1F(X( I) .AMAXX) 

00001940 

AMAXY  »  MAX1F(Y( I). AMAXY) 

00001950 

AMINX  *  MIN1F(X(I) .AMINX) 

00001960 

1020 

AMINY  =  MINIF(YU).  AMINY) 

00001970 

AMAXA  »  MAX1F(ABSF( AMAXX). 

ABSFJ AMAXY).  ABSF(AMINX).  ABSF( AMINY))  00001980 

IF(AMAXA  -  1.0E+99)1022. 1022. 1021                                   00001990 

1021 

PRINT  1102 

00002000 

1102 

FORMAT  (/»  58H  NO  X  OR  Y  VALUE  MAY  EXCEED  1.0E+99  IN  ABSOLUTE  MAGN00002010 

LITUDE.  ) 

00002020 

GO  TO  1005 

00002030 

1022 

IF(ABSF(AMAXX  -  AMINX)  -  1, 

.0E-97J1023.1025.1025                     00002040 

1023 

IF(ABSF(AMAXY  -  AMINY)  -  1, 

.0E-97) 1024,1025.1025                     00002050 

1024 

PRINT  1103 

00002060 

1103 

FORMAT  (/,  38H  ALL  POINTS  HAVE  THE  SAME  COORDINATES.  )             00002070 

GO  TO  1005 

00002080 

40 


1025 

IF( ITESTJ4.7.4 

4 

IF(MODCURV  -  2)5.240.5 

5 

IF(MODCURV  -  3)6.240t6 

6 

PRINT  101 

101 

FORMAT  </.  17H  ILLEGAL  MODCI 

GO  TO  1005 

7 

IF(MODCURV)6»9.8 

8 

IF(MODCURV  -  1)6*9*6 

9 

IF(IWIDE)10,11,12 

10 

ITIT  »  5HIWIOE 

PRINT  102»  ITIT.  ITIT 

102 

FORMAT  (/.   9H  ILLEGAL  .A5»i 

] 

I             5H  »  8.  ./) 

11 

JWIDE  «  8 

GO  TO  14 

12 

IF( IWIOE  -  9)13.13.10 

13 

JWIDE  ■  I WIDE 

14 

IFUHIGHU5.16.17 

15 

ITIT  -  5HIHIGH 

PRINT  102.  ITIT.  ITIT 

16 

JHIGH  «  8 

GO  TO  19 

17 

IFdHIGH  -  15)18.18.15 

18 

JHIGH  »  IHIGH 

19 

NODEXAX  «  MODEXAX 

IF(MODEXAX)20.27.21 

20 

ITIT"  8HMODEXAX. 

PRINT  104.  ITIT.  IX 

104 

FORMAT  (/.   9H  ILLEGAL  »A8. 

1         Al.   7HAX  =  0.  »/) 

NODEXAX  =  0 

GO  TO  27 

21 

IFtMODEXAX  -  2)27.22.20 

22 

IFdXUP  -  JHIGHJ24. 24.23 

23 

ITIT  =  8HIXUP. 

PRINT  104.  ITIT.  IX 

NODEXAX  »  0 

GO  TO  27 

24 

IF(IXUP)23.26.26 

26 

JXUP  •»  IXUP 

27 

NODEYAX  =  MODEYAX 

IF(MODEYAX)28»35.29 

28 

ITIT«8HMODEYAX. 

PRINT  104.  ITIT.  IY 

NODEYAX  *  0 

GO  TO  35 

29 

IF (MODEYAX  -  2)35.30.28 

30 

IFUYRIGHT  -  JWIDE)32.32.31 

31 

ITIT  ■  8HIYRIGHT. 

PRINT  104.  ITIT.  IY 

NODEYAX  «  0 

60  TO  35 

32 

IFUYRIGHT)  31.34.34 

GRAPH  WILL  BE  PLOTTED  WITH  .A5. 


32H  GRAPH  WILL  BE  PLOTTED  WITH  MODE. 


00002090 

00002100 

00002110 

00002120 

00002130 

00002140 

00002150 

00002160 

00002170 

00002180 

00002190 

00002200 

00002210 

00002220 

00002230 

00002240 

00002250 

00002260 

00002270 

00002280 

00002290 

00002300 

00002310 

00002320 

00002330 

00002340 

00002350 

00002360 

00002370 

00002380 

00002390 

00002400 

00002410 

00002420 

00002430 

00002440 

00002450 

00002460 

00002470 

00002480 

00002490 

00002500 

00002510 

00002520 

00002530 

00002540 

00002550 

00002560 

00002570 

00002580 

00002590 

00002600 

00002610 
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34 


WITH  YDATA. 


35 
2235 


50 


52 
53 

114 


56 
57 
58 


1030 

1031 
1104 


1032 
1033 
1039 
1034 

1036 


JYRIGHT  *  IYRIGHT 

INITIALIZE  PRIOR  TO  SCALING  AND  AXIS  LOCATING. 

IFLAG  =  0  FOR  PASS  WITH  XDATA.  IFLAG  »  1  FOR  PASS 
DO  2235  IOTA=l,12 
JJTITLE( IOTA)  =  ITITLE(IOTA) 
IFLAG  «  0 
BETA  »  0. 
SCALE  «  EXSCALE 
IAXIS  =  JYRIGHT 
MODE  *  NODEYAX 
ISI2E  *  JWIDE 
IXY  *  IX 
IYX  »  IY 
AMAX  =  AMAXX 
AMIN  ■  AMINX 
GO  TO  52 
IFLAG  »  1 
BETA  =  0. 
SCALE  =  YSCALE 
IAXIS  *  JXUP 
MODE   =  NODEXAX 
ISIZE  =  JHIGH 
AMAX  »  AMAXY 
AMIN  =  AMINY 
IXY  ■  IY 
IYX   »  IX 

CHECK  SCALE  AND  GO  TO  FIXED  OR  AUTO  SCALE  ROUTINES, 
IF(SCALE)53»59.56 
PRINT  114.  IXY»  IXY 

FORMAT  (/.   9H  ILLEGAL  •  Al ,39HSCALE. 
1TO  .Al,  7H-SCALE.  ,/) 
GO  TO  59 

EXPRESS  FIXED  SCALE  IN  E  FORMAT  WITH  ONE  FIGURE  SIGNIFICANCE, 
IF(SCALE  -  1.0E+99)57.53,53 
IF(SCALE  -  1.0E-99)53.53,58 
CALL  SCALE  IT (SCALE. I SCAL10. FACTOR. 1) 
SCALE  =  FACTOR*10.**ISCAL10 

CHECK  AND  COMPUTE  AXIS  LOCATION  IF  NECESSARY.  FIXED  SCALE 

CASE.  ITAG  =  0  IF  ORIGIN  ON  GRAPH  OR  1  IF  IT  IS  SUPPRESSED. 
IF(MODE  -  1)1032.1031.1030 
ITAG  =0 
GO  TO  203 

PRINT  1104  .  IYX.  IXY.  IXY 

FORMAT  (/.  5H  MODE . Al .24HAX  MUST  NOT  BE  1 
1  (AUTO-SCALE).  GRAPH  WILL  BE  PLOTTED  WITH 
GO  TO  59 

IF(ABSF(AMAX  -  AMIN)  -  1 .OE-97 ) 1033 . 1038 » 1038 
IF(ABSF(AMAX)  -  1 .OE-97 ) 1034 » 1039  ♦  1039 
IF (AMAX) 1036. 1034, 1037 
IAXIS  =  ISIZE/2 
GO  TO  1030 
IAXIS  =  ISIZE 
GO  TO  1030 


GRAPH  WILL  BE  PLOTTED  WITH 


UNLESS  .A1.57HSCALE  I 
AUTO  .Al,  7H-SCALE.  ♦ 


00002620 
00002630 
00002640 
00002650 
00002660 
00002670 
00002680 
00002690 
00002700 
00002710 
00002720 
00002730 
00002  740 
00002750 
00002760 
00002770 
00002780 
00002790 
00002800 
00002810 
00002820 
00002830 
00002840 
00002850 
00002860 
00002870 
00002880 
00002890 
00002900 
AUOO002910 
00002920 
00002930 
00002940 
00002950 
0000296( 
0000297( 
0000298C 
0000299( 
0000  3  00( 
000030K 
0O00302( 
0000303( 
00003040 
S  000003050 
/)  00003  060 
00003070 
00003080 
00003090 
00003100 
00003110 
00003120 
00003130 
00003140 
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1037  IAXIS  =  0  00003150 
GO  TO  1030  00003160 

1038  IF(SIGNF( 1..AMAX)  -  S IGNF ( 1 • .AMIN )) 1040 > 1039 ♦ 1040  00003170 
1040  ASIZE  =  ISIZE  00003180 

IAXIS  =  -AMIN/(AMAX  -  AMIN)*ASIZE  +0,5  00003190 

GO  TO  1030  00003200 

AUTO  SCALE  ROUTINE.  00003210 

59  IF (MODE  -  1)60.64.69  OO003220 

60  AMAX  =  MAXIFtO.t  AMAX)  00003230 
AMIN  =  MINlF(0.t  AMIN)  00003240 

64  IF(ABSF(AMAX  -  AMIN)  -  1 .0E-97 ) 65 #68 .68  00003250 

65  PRINT  116t  IXY»  IXY.  IYX  00003260 
116  FORMAT  (/.   5H  ALL  .A1.47H  VALUES  EQUAL.  AUTO  SCALE  POSSIBLE  ONLY  00003270 

1IF  THE  »A1»29H  VALUES  ARE  NON-ZERO  AND  MODE. Alt  7HAX  =  2.  )  00003280 

GO  TO  1005  00003290 

68  ASIZE  =  ISIZE  00003300 
SCALE  =  (AMAX  -  AMIN) /ASIZE  00003310 
GO  TO  83  00003320 

69  IF(ABSF(AMAX  -  AMIN)  -  1 .0E-97 ) 70.74.74  00003330 

70  IF(ABSF(AMAX)  -  1 .0E-97 ) 71 .74. 74  00003340 

71  PRINT  118.  IXY  00003350 
118  FORMAT  (/.   5H  ALL  »A1»38H  VALUES  ZERO.  AUTO  SCALE  NOT  POSSIBLE.  J00003360 

GO  TO  1005  00003370 

74  IF(AMAX  -  1.0E-97)  76»75.75  00003380 

75  IFUSIZE  -  IAXISJ77.76.77  00003390 

76  SCALE1  =  0.  00003400 
GO  TO  78  00003410 

77  AXIS  =  IAXIS  00003420 
ASIZE  ■  ISIZE  00003430 
SCALE1  =  AMAX/( ASIZE  -  AXIS)  00003440 

78  IF(AMIN  +  1.0E-97J79.79.30  00003450 

79  IF(  IAXISJ81.80.81  00003460 

80  SCALE2  =  0.  00003470 
GO  TO  82  00003480 

81  AXIS  =  IAXIS  00003490 
SCALE2  =  -AMIN/AXIS  00003500 

82  IF(SCALE1  +  SCALE2 ) 1984, 1982 . 1984  00003510 

1982  PRINT  1983.  IYX.  IYX  00003520 

1983  FORMAT  (/»  56H  NONE  OF  THE  PLOT  LIES  ON  THE  GRAPH  WITH  THIS  SPECIF00003530 
1IED  .A1.47H-AXIS  LOCATION.  GRAPH  WILL  BE  PLOTTED  WITH  M0DE.A1.  00003540 
2   7HAX  =0.   ./)  00003550 

MODE  =  0  00003560 

GO  TO  60  00003570 

1984  SCALE  =  MAX1F ( SCALE1 .  SCALE2)  00003580 

83  CALL  SCALEIT(SCALE.  ISCAL10. FACTOR .3 )  00003590 
IF(FACT0R  -  5.05)85.85.84  00003600 

84  FACTOR  =  1  00003610 
ISCAL10  ■  ISCAL10  +  1  00003620 
GO  TO  90  00003630 

85  IFIFACTOR  -  2.02)87.87.36  00003640 

86  FACTOR  =  5  00003650 
GO  TO  90  00003660 

87  IF(FACTOR  -  1.01)89.89.88  00003670 
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88 

89 
90 


91 
92 
93 

94 
95 


96 
120 


1 


99 

97 
200 

201 

202 


203 
204 


205 


FACTOR  =  2 

GO  TO  90 

FACTOR  =  1 

SCALE  =  FACTOR«10.**ISCAL10 

COMPUTE  AXIS  LOCATION  IF  NECESSARY,  AUTO  SCALE  CASE, 
IF (MODE  -  1)92.91*93 

IF(SIGNF( 1..AMAX)  -  S TGNF ( 1 • » AMIN ) ) 92.94.92 
IAXIS  «  -AMIN/SCALE  +  0.5 
I  TAG  *  0 
GO  TO  203 
IF(AMAX)95»95.200 
IAXIS  »  ISIZE 
BETA  *  -AMAX/SCALE 
IF(BETA  -  l.E+12)99. 99.96 
PRINT  120.  IXY 

FORMAT  (/.  15H  THE  ORIGIN  OF  ,A1.  43H  CANNOT  BE  OFFSET  MORE  THAN 
•OE+12  INCHES.  ) 
GO  TO  1005 
I8ETA  -  BETA  ♦  0.5 
BETA  =  -IBETA 

BETA  IS  THE  NUMBER  OF  INCHES  OF  ORIGIN  SUPPRESSION. POS I T IVE  I 

TRUE  ORIGIN  IS  BELOW  OR  TO  LEFT  OF  THE  GRAPH. 
IFfBETA  +  D97.97.93 
ITAG  =  1 
GO  TO  203 
IAXIS  «  0 
BETA  =  AMIN/SCALE 
IF(BETA  -  l.E+12)201. 201.96 
IBETA  =  BETA  +  0.5 
BETA  *  IBETA 
IF(BETA  -  l.)93,202.202 
ITAG  =  1 

RELEASE  RESULTS  TO  REMAINING  PART  OF  PROGRAM.  START  SECOND 

PASS  FOR  Y  VALUES  IF  NOT  YET  COMPUTED. 
IF( IFLAG)205»204.205 
SCALEX  =  SCALE 
IXFACT   =  FACTOR 
IXSC10  =  ISCAL10 
IXAXIS  =  IAXIS 
ITAGX   =  ITAG 
ISIZEX  =  ISIZE 
BETAX  »  BETA 


C 

c 


GO  TO  50 
BETAY  »  BETA 
SCALEY  «  SCALE 
IYFACT  =   FACTOR 
IYSC10  =  ISCAL10 
IYAXIS  =  IAXIS 

NOW  WRITE  RECORDS. 
ITAGY   »  ITAG 
ISIZEY  »  ISIZE 

THIS  COMPLETES  CALCULATION  OF  SCALE  FACTORS  ETC. 

TAPE  RECORDS.  FIRST.  THE  SCALE  FACTOR  TITLES. 


NOW  GENERATE 


00003680 
00003690 
00003700 
00003710 
00003720 
00003730 
00003740 
00003750 
00003760 
00003770 
00003780 
00003790 
00003800 
00003810 
00003820 

100003830 
00003840 
00003850 
00003860 
00003870 

F00003880 
00003890 
00003900 
00003910 
00003920 
00003930 
00003940 
00003950 
00003960 
00003970 
00003980 
00003990 
00004000 
00004010 
00004020 
00004030 
00004040 
00004050 
00004060 
00004070 
00004080 
00004090 
00004100 
00004110 
00004120 
00004130 
00004140 
00004150 
00004160 
00004170 
00004180 
00004190 
00004200 
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206  JXTIT( 1)  =  8H1      X- 
JXTIT(2)  =  8HSCALE  = 
JXTITO)  =  ICODE(SCALEX) 
JXTIT(4)  =  8H  UNITS/I 
JXTIT(5)  ■  8HNCH. 
JYTIT(l)  *  8H1      Y- 
JYTIK2)  ■  8HSCALE  * 
JYTITO)  «  ICODE(SCALEY) 
JYTIT(4)  «  8H  UNITS/I 
JYTIT<5)  *  8HNCH. 
DO  9206  1=6.11 
JXTITt I)  «  8H 
JYTIT( I)  »  8H 
IF(ITAGX)211.211.207 
IF(BETAX)208. 208*209 


9206 

207 
208 

209 
210 


JXTIT(7) 

GO  TO  210 

JXTIT(7) 

JXTIK8) 

JXTIK9) 

JXTIT(IO) 

JXTIT(ll) 


8H 


ADO  - 


211 
212 
213 

214 
215 


216 


8H    ADD  + 

ICODE  (BETAY*SCALEY) 

8H  UNITS  T 


8H    ADD  ♦ 

ICODE  (BETAX*SCALEX) 

8H  UNITS  T 

8HO  ALL  X 

8HVALUES. 
IF( ITAGY>216.216.212 
IF(BETAY)213»213.214 
JYTIT(7)  «  8H    ADD  - 
GO  TO  215 
JYTIT(7) 
JYTIT(8) 
JYTIT(9) 

JYTIT< 10)*8HO  ALL  Y 
JYTIT(H)-  8HVALUES, 
ICONT(l)  -  ICONTRL  +  4 
C  INSERT  TITLE  SIZE  (02B)  AHFAD  OF  MAIN  TITLE  RECORD. 

CALL  ISHIFT6  (ITITLE.  LTITLE) 
C         TEST  FOR  ALL  BLANK  TITLES. 
ICHECK  =  8H 
DO  9075  1=1.6 

IF( ITITLE( I )-ICHECK)  9074»9075 .9074 
IF(ITITLE(U  )  9080.9075.9080 
CONTINUE 
IT1  »  1 

ICONT(l)  -  ICONT(l)  -  1 
GO  TO  9081 
IT1  =  0 

DO  9085  1*7.12 

IF  (ITITLE(I)  -  ICHECK)  9084.9085.9084 
IF  (  ITITLEt I) )9090, 9085. 9090 
CONTINUE 
IT2  =  1 

ICONT(l)  »  ICONT(l)  -  1 
GO  TO  9091 
IT2  =  0 

NOW  GENERATE  AXES  RECORDS. 


9074 
9075 


9080 
9081 

9084 
9085 


9090 


00004210 
00004220 
00004230 
00004240 
00004250 
00004260 
00004270 
00004280 
00004290 
00004300 
00004310 
00004320 
00004330 
00004340 
00004350 
00004360 
00004370 
00004380 
00004390 
00004400 
00004410 
00004420 
00004430 
00004440 
00004450 
00004460 
00004470 
00004480 
00004490 
00004500 
00004510 
00004520 
00004530 
00004540 
00004550 
00004560 
00004570 
00004580 
00004590 
00004600 
00004610 
00004620 
00004630 
00004640 
00004650 
00004660 
00004670 
00004680 
00004690 
00004700 
00004710 
00004720 
00004730 
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9091 


240 
9007 
9700 

9701 
242 


9241 
9242 


LFTMGN  =  0*100 

IBOTMGN  *    1*100 

IH  =  LFTMGN 

JH  =  IBOTMGN  +IYAXIS*100 

LH  =  ISIZEX*100 

IHL  *  LFTMGN  +  ISIZEX*100  -  107 

KAXIS(l)  =  lPACK12(IH»JH.LHf IHL) 

JHL  =  JH  -  13 

IHL2  =  -100 

IVH   =  (ISIZEX  -IXAXIS  -  1)*IXFACT 

IVH2  =  -IXFACT 

KAXIS(2)  =  IPACK12(JHL.IHL2.IVH*IVH2> 

NH  =  ISIZEX 

ISH  =  8H       14 

IV  =  LFTMGN  +  IXAXIS*100 

JV  =  IBOTMGN 

KAXISO)  =  IPACK12(NH»ISH»IV»JV) 

LV  =  ISIZEY*100 

IVL  =  IV  -  3 

JVL  «  IBOTMGN  +  ISIZEY*100  -  107 

JVL2  *  -100 

KAXIS(4)  =  IPACK12(LV*IVL»JVL,JVL2) 

IVV  =  (ISIZEY  -  IYAXIS  -  1)*IYFACT 

IVV2  =  -IYFACT 

INV  =  ISIZEY 

ISV  =  8H       11 

KAXIS<5)  «  IPACK12(IVV»IVV2»INV»ISV) 

NOW  GENERATE  CURVES. 
SCX  =  100./SCALEX 
SCY  =  100./SCALEY 
EXAXIS=  IXAXIS*100 
YAXIS  *    IYAXIS*100 
ALFTMGN  =  LFTMGN 
BOTMGN  =  IBOTMGN 

SHIFTX  *EXAXIS  -  BETAX*100.  +  ALFTMGN 
SHIFTY  =  YAXIS  -  BETAY*100,  +  BOTMGN 
EXSIZE=  ISIZEX*100  +  LFTMGN  +  60 
SIZEX  *  LFTMGN  -  60 
YSIZE  =  ISIZEY*100  +  IBOTMGN  +  70 
SIZEY  »  IBOTMGN  -  70 
ICURV(l)  =  0 

IF( IPOINT)9010. 9007*9010 
IF(XMODF(NUMPTS»2) ) 9700 t9701 #9700 
I  SWITCH  =  1 
GO  TO  242 
I  SWITCH  *  2 
INUM  «  (NUMPTS  +  l)/2 
DO  244  1=1* INUM 
CI  =  X(2*I-1)*SCX  +  SHIFTX 
C2  =  Y(2*I-1)*SCY  +  SHIFTY 
IF( I-INUM)241»9241»241 
GO  TO  (9242.241)#ISWITCH 
C3  =  CI 


00004740 
00004750 
00004760 
00004770 
00004780 
00004790 
00004800 
00004810 
00004820 
00004830 
00004840 
00004850 
00004860 
00004870 
00004880 
00004890 
00004900 
00004910 
00004920 
00004930 
00004940 
00004950 
00004960 
00004970 
00004980 
00004990 
00005000 
00005010 
00005020 
00005030 
00005040 
00005050 
00005060 
00005070 
00005080 
00005090 
00005100 
00005110 
00005120 
00005130 
00005140 
00005150 
00005160 
00005170 
00005180 
00005190 
00005200 
00005210 
00005220 
00005230 
00005240 
00005250 
00005260 
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241 
9243 


244 
246 

9010 
247 

1260 

260 

261 

265 
9268 

9269 

9270 

9271 

9015 
270 

9020 
272 

273 


1273 


1274 
1275 


C4 
GO 
C3  ■ 
C4  ' 
CI  « 
IC1' 
C2  ■ 
IC2  = 
C3  i 
IC3  = 
C4  ' 
IC4' 


IDUMMY) 


1) 


=  C2 

TO  9243 

=  X(2*I)*SCX  +  SHIFTX 

«  Y(2»I )*SCY  +  SHIFTY 

=  MINIF(CI.EXSIZE) 

MAX1F(C1.  SIZEX) 

MIN1F(C2.  YSIZE) 

MAX1F(C2»  SIZEY) 

MIN1F(C3.EXSIZE) 

MAX1F(C3»  SIZEX) 

MIN1F(C4.  YSIZE) 

MAX1FCC4.  SIZEY) 
ICURVU+1)  «  IPACK12< IC1.IC2.IC3.IC4) 
II  =  INUM  ♦  3 

CALL  IPACKLKLABEL.  LABEL1» 
ICURV( 1 1  — 1 )  *  LABEL1 
ICURV(II)  *  ICURV4 
IF(MODCURV  -  l)247t247.9015 
CALL  IREADY  (IDUMMY) 
IF( IDUMMY)5000.1260.5000 
CALL  IWRITE  (ICONT.  IDUMMY. 
IF( I DUMMY) 5000. 260 • 5000 
CALL  IWRITE  (JXTIT.  IDUMMY.ll) 
IF( IDUMMYJ5000. 261. 5000 
CALL  IWRITE  (JYTIT.  IDUMMY.ll) 
IF( IDUMMYJ5000. 265. 5000 
IF(IT1)9269, 9268. 9269 
CALL  IWRITE  (LTITLE.  IDUMMY.  7) 
IF( IDUMMY)5000. 9269, 5000 
IF( IT2)9271»9270»9271 
CALL  IWRITE  (LTITLE(7).  IDUMMY. 
IF( IDUMMY)5000.9271.5000 
CALL  IWRITE  (KAXIS.  IDUMMY. 
IF( IDUMMY)5000. 9015. 5000 
IF( IPOINTJ9020. 270. 9020 
CALL  IWRITE  (ICURV.  IDUMMY. 
IF( IDUMMY)50O0. 9020. 5000 
IFfMODCURV  -  1)272.272.9025 
IFUGRID  -  1)9025.273.9025 

GENERATE  GRID  IF  CALLED  FOR. 
1X100  =  ISIZEX*100 
IY100  »  ISIZEY*100 
NEXT1  «  IBOTMGN 
NEXT2  »  LFTMGN  +  1X100 
JGRID(l)  *  0 
DO  1274  J=1.19.2 
JGRIDU+1)  ■  IPACK12 
IF(NEXT1  -  IBOTMGN  - 
NEXT1  =  NEXT1  +  100 
JGRIDU+2)  =  IPACK12 
IF(NEXT1  -  IBOTMGN  - 
NEXT1  =  NEXT1  +  100 
JGRID(J+2)  =  IPACK12 


7) 


5) 


II) 


(LFTMGN.  NEXT1.  NEXT2.  NEXT1) 
IY10O1273.  1275.  1275 


(NEXT2.  NEXT1.  LFTMGN. 
IY100)1274. 1276. 1276 


NEXT1) 


(NEXT2.  NEXT1.  NEXT2.  NEXTl) 


00005270 
00005280 
00005290 
00005300 
00005310 
00005320 
00005330 
00005340 
00005350 
00005360 
00005370 
00005380 
00005390 
0O005400 
00005410 
00005420 
00005430 
00005440 
00005450 
00005460 
00005470 
00005480 
00005490 
00005  500 
00005510 
00005520 
00005530 
00005540 
00005550 
00005560 
00005570 
00005580 
00005590 
00005600 
00005610 
00005620 
00005630 
00005640 
00005650 
00005660 
00005670 
OTO05680 
00005690 
00005700 
00005710 
00005720 
00005730 
00005740 
00005750 
00005760 
00005770 
00005780 
00005790 


hi 


1276 


1277 


1278 


1279 
1280 
1281 


9025 


9030 


9031 
9032 
9033 
9034 

9035 


9036 


9037 


9038 


9039 


JGRIO(J+3)  =  ICURV3 

JGRIDU+4)  =  ICURV4 

CALL  IWRITE( JGRID.  IDUMMY,  J+4) 

IF( IDUMMY)5000,1277,5000 

NEXT1  =  LFTMGN 

NEXT2  =  IBOTMGN  +  IY100 

DO  1279  J=l,ll»2 

JGRIDU+1)  =  IPACK12  (NEXTl,  IBOTMGN.  NEXT 

IF(NEXT1  -  LFTMGN  -  1X100)1278.1280,1280 

NEXTl  =  NEXTl  +  100 

(NEXTl,  NEXT2.  NEXTl. 
1X100)1279.1281,1281 


(NEXTl,  NEXT2,  NEXTl, 


J+4) 


JGRIDU+2)  =  IPACK12 
IF(NEXT1  -  LFTMGN  - 
NEXTl  =  NEXTl  +  100 
JGRID(J+2)  =  IPACK12 
JGRID(J+3)  =  ICURV3 
JGRID(J+4)  =  ICURV4 
CALL  IWRITE  (JGRID,  IDUMMY, 
IF( IDUMMY)5000, 9025,5000 
IF (I  POINT) 9030 ,276, 9030 

GENERATE  POINT  PLOT  RECORDS 
IOUT  =  0 

CALL  IPACKL1  (LABEL,  LABEL1,  LABEL2) 
DO  9050  I=1,NUMPTS 
CI  =  X(  I)*SCX  +  SMIFTX 
C2  =  Y( I )»SCY  +  SHIFTY 
IF(C1  -  EXSIZE)9031, 9031, 9034 
IF(C2  -  YSIZE)9032, 9032, 9034 
IF(C1  -  SIZEX)9034, 9033,9033 
IF(C2  -  SIZEY)9034, 9035, 9035 
IOUT  =  IOUT  +1 
GO  TO  9050 
IC1  =  CI 
IC2  =  C2 
GO  TO  (9036,9037.9038,9039,9040) .IPOINT 

GENERATE  CROSS, 
ICURV(2)  =  IPACK12  (IC1-5, 

(  IC1   , 
( IC1+5, 


1,  NEXT2) 
IBOTMGN) 
NEXT2) 


IF  CALLED  FOR, 


IPACK12 
IPACK12 


IC2-5, 
IC2  . 
IC2-5. 


ICURV(3) 
ICURV(4) 
GO  TO  9041 

GENERATE  PLUS. 
ICURV(2)  ■  IPACK12  (  IC1   • 
ICURV(3)  =  IPACK12  ( IC1   » 
ICURV(4)  =  IPACK12  (IC1+5. 
GO  TO  9041 

GENERATE  SQUARE. 
ICURV(2)  =  IPACK12  (IC1+4, 
ICURV(3)  =  IPACK12  (IC1-4. 
ICURV(4)  =  IPACK12  (IC1+4. 
GO  TO  9041 

GENERATE  DIAMOND. 
ICURV(2)  =  IPACK12  (IC1+5.  IC2 
ICURV(3)  «  IPACK12  (IC1-5.  IC2 
ICURV(4)  *  IPACK12  (ICl-^5.  IC2 


IC1+5. 
IC1-5. 

rci+5. 


IC2-5, 
IC2  ♦ 
IC2   ♦ 


IC2-4. 
IC2+4. 
IC2-4. 


IC1  , 
IC1-5, 
IC1+5, 


IC1+4, 
IC1-4, 
IC1+4, 


IC1  , 
IC1  , 
IC1+5, 


C2  +  5) 
C2  +  5) 
C2-5) 


C2  +  5) 
C2  ) 
C2   ) 


C2+4) 
C2-4) 
C2-4) 


C2  +  5) 
C2-5) 
C2   ) 
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00006060 
00006070 
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9040 


9041 
9042 

9043 
9044 

9045 
9046 


9050 

9048 
9104 

276 
277 
278 

130 


656 
279 


5000 
5001 

5002 
5100 


5670 


1NEX 


GO  TO  9041 

GENERATE  TRIANGLE. 
ICURV(2)  =  IPACK12  (IC1+5.  IC2-3t  IC1   .  IC2+6) 
ICURV(3)  =  IPACK12  (IC1-5.  IC2-3.  IC1+5.  IC2-3) 
ICURV14)  ■  ICURV(3) 
IF ( I  -  NUMPTS)9043»9042.9043 
ICURV(5)  =  LABEL2 
GO  TO  9046 

IF( I  -  1)9045»9044»9045 
ICURV(5)  *  LABEL1 
GO  TO  9046 
ICURVC5)  =  ICURV3 
ICURV(6)  =  ICURV4 
CALL  IWRITE  (ICURV.  IDUMMY.  6) 
IF( IDUMMY)5000.9050»5000 
CONTINUE 

IF (IOUT) 9048 t276» 9048 
PRINT  9104,  IOUT 
FORMAT  (/.  1X»  12.  29H  POINT(S)  WERE  OFF  THE  GRAPH.  ,/) 

SET  UP  RETURN. 
IF (MODCURV) 277.278 ,277 
IF(MO0CURV  -  3)279.278.279 
ITEST  *  0 

PRINT  130,( JJTITLE(I) .1=1.12) 

FORMAT  (/.  19H  GRAPH  TITLED  .  .   . 6A8 »/ , 19X . 6A8 , 
L  24H   •  •  HAS  BEEN  PLOTTED.  ./.1H0) 

IDUMMY  =  ITYP2( IDUMMY) 
IF( IDUMMYJ5670. 656. 5670 
LAST  =  0 
RETURN 
ITEST  «  1 

IDUMMY  »  I CLOCK (IDUMMY) 
LAST  =  0 
RETURN 

THESE  ARE  THE  NORMAL  RETURNS. 

NOW  SET  UP  THE  RETURN  FOLLOWING  A  TAPE  ERROR. 
IF(MODCURV  -  1)5001.5001.5002 
IDUMMY  *  ITYPEK  IDUMMY) 
GO  TO  247 
PRINT  5100 

FORMAT  (/.  36H  TAPE  ERROR  IN  WRITING  GRAPH  OUTPUT.  ) 
IDUMMY  ■  ITYPEK IDUMMY) 
GO  TO  1007 

IDUMMY  «  ITYPEK  IDUMMY) 
END 


SUBROUTINE  I  READY  (IDUMMY) 

SELECTS  TAPE  8  (WILL  LOOP  UNTIL  READY).  WRITES  EOF  ON  8. 
MACHINE  LANGUAGE  WILL  NOT  BE  NECESSARY  IN  FORTRAN  62-3. 
LOCdFIVE  *  5). 
EXF  (52041B)     EXF7  (52000B).    SELECT  READ  AND  WAIT  TAPE. 
EXF7(00050B)     SLJ  (1RDY).       EXIT  ON  CH  5  ACTIVE. 
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EXF7(52000B) 

SLJ  (1NEX). 

EXIT  ON  TAPE  READY. 

00006860 

LDA  ( IFIVE) 

SAU  (1BUF). 

TERMINATE 

00006870 

1BUF 

EXF5(N). 

BUFFER. 

00006880 

1RDY 

EXF  (52041B) 

EXF7(52000B). 

SELECT  AND  WAIT  TAPE  8. 

00006890 

ENA  (0), 

CLEAR  A. 

00006900 

EXF  (02000B) 

EXF  (52006B). 

STOP  CLOCK  AND  BACKSPACE  8. 

00006910 

EXF7(52001B) 

SLJ  (1END). 

EXIT  IF  NOT  AT  LOAD  POINT. 

00006920 

-EXF7(52000B). 

WAIT  TAPE  8. 

00006930 

EXF7(52007B) 

SLJ  (1E0F). 

EXIT  IF  NO  EOF. 

00006940 

ENA  (I DUMMY) 

SAU  (2BUF). 

MOVE 

00006950 

INA  (1) 

STA  (IFIVE). 

FORWARD 

00006960 

2BUF 

EXF5(N) 

EXF7(52000B). 

OVER  RECORD. 

00006970 

1E0F 

ENA  (0) 

EXF7(00061B). 

CLEAR  A,  WAIT  CH  6. 

00006980 

EXF  (62041B) 

EXF7(62000B). 

SELECT  WRITE  AND  WAIT  TAPE. 

00006990 

EXF  (62041B) 

EXF7(62000B). 

SELECT  AND  WAIT  TAPE  8. 

00007000 

EXF  (62003B) 

EXF7(62000B). 

WRITE  EOF  AND  WAIT. 

00007010 

EXF7(62007B) 

ENA(IO). 

EXIT  ON  NO  END  OF  TAPE. 

00007020 

1END 

STA  (IDUMMY). 
END 

00007030 
00007040 

C 

00007050 

SUBROUTINE  IWRITE ( ISTART,  IDUMMYi 

»  IWORDS) 

00007060 

C 

WRITE  RECORD  OF  IWORDS»  STARTING  WITH  ISTART.  PUT  IDUMMY  «  0 

00007070 

C 

IF  RECORD 

CORRECTLY  WRITTEN* 

OTHERWISE  SET  NON-ZERO. 

00007080 

* 

MACHINE  LANGUAGE  WILL  NOT  BE 

NECESSARY  IN  FORTRAN  62-3. 

00007090 

LOCdSIX  =  6). 

00007100 

-EXF7(0O061B). 

WAIT  CH  6. 

00007110 

EXF  (62041B) 

EXF7(62000B). 

SELECT  WRITE»  WAIT  TAPE. 

00007120 

EXF  (62041B) 

EXF7(62000B). 

SELECT  AND  WAIT  TAPE  8. 

00007130 

ENQ  (111B). 

SET  COUNTER. 

00007140 

1AGN 

ENA  (ISTART) 

INA(l). 

STARTING  ADDRESS. 

00007150 

SAL  (1BUF) 

ADD( IWORDS). 

TERMINAL  ADDRESS. 

00007160 

1BUF 

STA(ISIX) 

EXF6(N). 

BUFFER  OUT. 

00007170 

ENA(O) 

EXF7(62000B). 

CLEAR  A.  WAIT  TAPE  8. 

00007180 

EXF7(62007B) 

SLJ  (1END). 

EXIT  IF  NO  END  OF  TAPE. 

00007190 

EXF7(62003B) 

SLJ  (2AGN). 

EXIT  IF  NO  PARITY  ERROR. 

00007200 

EXF7(62004B) 

SLJ  (2END). 

EXIT  IF  LENGTH  ERROR. 

00007210 

2AGN 

EXF  (62006B) 

EXF7(62000B). 

BACKSPACE  AND  WAIT. 

00007220 

ORS  (3) 

QJPK1AGN). 

TRY  WRITE  3  TIMES. 

00007230 

1END 

EXF  (62003B) 

ENA  (10). 

WRITE  EOF,  NON  ZERO  A. 

00007240 

2END 

STA  (IDUMMY). 
END 

STORE  RESPONSE. 

00007250 
00007260 

C 

00007270 

FUNCTION  ITYP2 

( IDUMMY) 

00007280 

C 

TYPE  WORD 

GRAPH. 

00007290 

* 

WILL  NEED 

REWRITING  IN  FORTRAN  62-3. 

00007300 

CON(LC  *  57B, 

Ml  =  4513123015050000B  )• 

00007310 

LOCUTWO  ■  2). 

00007320 

-EXF7(00061B). 

WAIT  CH  6. 

00007330 

EXF  (62041B) 

EXF7(62000B). 

SELECT  AND  WAIT  TAPE. 

00007340 

EXF  (62041B) 

EXF7(62000B). 

SELECT  AND  WAIT  TAPE  8. 

00007350 

EXF  (62003B) 

EXF  (01000B). 

WRITE  EOF.  START  CLOCK. 

00007360 

ENA  (0) 

EXF7(62000B). 

CLEAR  A.  WAIT  TAPE  8. 

00007370 

EXF7(62007B) 

ENA  (10). 

EXIT  IF  NO  END  OF  TAPE. 

00007380 
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STA  (ITYP2). 

-EXF7C00021B). 
EXF7(11141B) 
EXF  (21100B) 
STA  (ITWO) 

-EXF7(00021B). 
1TYP   EXF  (2110OB) 
STA  ( ITWO) 

END 


SLJ  (1TYP). 
ENA  (LC+1). 
EXF2(LC). 

ENA  (Ml+1). 
EXF2(M1). 


STORE 

•  RESPONSE. 

WAIT 

CH 

2. 

EXIT 

IF 

UPPER  C 

TYPE 

LOWER 

CASE. 

WAIT 

CH 

2. 

TYPE 

GRAPH 

CASE. 


1TYP 


FUNCTION  ITYPE1  (IDUMMY) 

REWIND  TAPE  8.  REQUEST  NEW  TAPE.  AND  WAIT  TILL  READY, 
WILL  NEED  REWRITING  IN  FORTRAN  62-3. 
CON(LC  *  57B.  Ml  =  451 5112030242004B .  M2 
1     M3  =  0401301520043342B). 
RSVCMESS  *  3). 
LOCUTWO  =2). 
-EXF7(00061B). 
EXF  (62041B) 
EXF  (62041B) 
EXF  (62003B) 
EXF  (62007B) 
EXF7(11141B) 
EXF  (21100B) 
( ITWO) 
(Ml) 
(M2) 
(M3) 
-EXF7(00021B). 
EXF  (211O0B) 
STA  (ITWO) 
-EXF7(00061B). 
-EXF7(62000B). 
EXF  (62041B)     EXF7 ( 62000B ) .     SELECT  AND  WAIT  TAPE  8 
EXF  (OlOOOB). 
END 


STA 
LDA 
LDA 
LDA 


EXF7(62000B) 
EXF7(62000B) 
EXF7(62000B) 
EXF7(00021B) 
SLJ  (1TYP). 
ENA  (LC+1). 
EXF2(LC). 
STA  (MESS). 
(MESS+1) 
(MESS+2) 


STA 
STA 


ENA  (MESS+3) 
EXF2JMESS). 


»  1103302204062031B, 


WAIT  CH  6. 

SELECT  AND  WAIT  TAPE. 

SELECT  AND  WAIT  TAPE  8. 

WRITE  EOF  AND  WAIT. 

REWIND  WITH  INTERLOCK*  WAIT  CH 

EXIT  IF  UPPER  CASE. 

TYPE 

LOWER  CASE. 


WAIT  CH  2. 
TYPE 

MESSAGE. 
WAIT  CH  6. 
WAIT  TAPE. 
SELECT  AND  WAIT 
START  CLOCK. 


FUNCTION  ICODE  (ANUMBER) 

CODES  ABSOLUTE  VALUE  OF  A  FLOATING  POINT  NUMBER  (BETWEEN 
1.0E-100  AND  1.0E+10O)  INTO  8-CHARACTER  BCD  WORD  OF  THE  FORM 
1.23E+45.  ICODE  =  8H0.00E+00  IF  MAGNITUDE  OUT  OF  RANGE. 
CHECK  AVAILABILITY  OF  LIBRARY  FUNCTIONS  IN  FORTRAN  62-3. 

DIMENSION  11(8) 

BNUMBER  *  ABSF( ANUMBER) 

IFJBNUMBER  -  1 .OE+lOO ) 7 »6 . 6 

IF(BNUMBER  -  1 .OE-lOO ) 6.6 . 2 

ICODE  =  8H0.00E+00 

RETURN 

THIS  IS  ERROR  EXIT. 

CALL  SCALEIT  (BNUMBER.  ISCAL10.  FACTOR.  3) 

ISIGNSC  =  XSIGNF(l.ISCALlO) 

ISCAL10  =  XABSF(ISCALIO) 

IFACT  «  FACTOR*100.001 
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11(8)  »  XMODF( ISCALIO, 10) 

II(7)*ISCAL10/10 

IF( ISIGNSC)4,3»3 


11(6) 
GO  TO 
11(6) 
11(5) 
11(4) 
11(3) 
11(2) 
11(11 
CALL 


«  8H        + 
5 

■  8H 

=  8H        E 

■  XMODF( IFACT,10) 

=»  (XMODF( IFACT,100) )/10 

■  8H        • 

■  IFACT/lOO 
IPACK  (II,  IPACKED) 


I CODE  » 

RETURN 

END 


IPACKED 


SUBROUTINE  SCALEIT  < ANUMBER, ISCALIO, FACTOR, MODE ) 
C  FINDS  FACTOR  (BETWEEN  1.0  AND  9.99...)  AND  SCALE  OF  10  AS 

C  DEFINED  BY      ANUMBER  =  FACTOR»10.*»ISCALlO. 

C  MODE  IS  THE  NUMBER  OF  SIGNIFICANT  FIGURES  REQUIRED.  THIS  MUST 

C  BE  BETWEEN  1  AND  10  OR  IT  WILL  BE  PUT  EQUAL  TO  SIX. 

»  CHECK  AVAILABILITY  OF  LOG10F  IN  FORTRAN  62-3. 

ISCAL10«LOG10F( ANUMBER) 

FACTOR  «  ANUMBER/10.**ISCAL10 

IF(FACTOR  -  0.1)1,2,2 

1  FACTOR  «  FACTOR«100. 
ISCALIO  *  ISCALIO  -  2 
GO  TO  8 

2  IF(FACTOR  -  1.0)3,8,4 

3  FACTOR  «  FACTOR*10. 
ISCALIO  *  ISCALIO  -  1 
GO  TO  8 

4  IF(FACTOR  -  100.0)6,5,5 

5  FACTOR  «  FACTOR/100. 
ISCALIO  *  ISCALIO  +  2 
GO  TO  8 

IF(FACTOR  -  10.0)8,7,7 
FACTOR  -  FACTOR/10. 
ISCALIO  =  ISCALIO  +  1 
IF(MODE)9,9,10 
MODE  =  6 
GO  TO  11 

IF (MODE  -  10)11,11,9 

IFACTOR  «  F ACTOR* 10. ♦* (MODE  -  1)  +  0.5 
FACTOR  «  IFACTOR 

FACTOR  *  FACTOR/10. **( MODE  -  1) 
IF(FACTOR  -  10.)13,12,12 
FACTOR  ■  1#  ■ 
ISCALIO  «  ISCALIO  +  1 
RETURN 
END 

SUBROUTINE  IPACK  (II,  IPACKED) 


6 

7 

8 
9 

10 
11 


12 


13 
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1NEX 
2NEX 


1NEX 


2NEX 


TAKES  8  SIX-BIT  WORDS  AND  PACKS  THEM  LEFT  TO  RIGHT 
IN  IPACKED.  IF  WORD  IS  ZERO.  12B  IS  SUBSTITUTED. 
CONVERT  TO  CODAP  FOR  FORTRAN  62-3, 
CONUZERO  *  12B). 


SIUK  ISAVE) 
LDAK  II) 
LDA  (I ZERO). 
LRS  (6) 
ISKK-1) 
STQ  ( IPACKED) 
END 


ENIK8). 
AJPK2NEX). 

INIK-2). 
SLJ  (1NEX). 
LIUK  ISAVE). 


SUBROUTINE  ISHIFT6  (IT 

INSERTS  02B  AHEAD 

WILL  HAVE  TO  BE  CON 

WATCH  ARRAY  INDEX  I 

CON< IBLANK  «  202020202 


SIUK  ISAVE) 
ENA  (2). 
LDQK  ITITLE) 
STAKLTITLE) 
ISKK6) 
LDQ  ( IBLANK) 
ENIK7) 
ENA  (2). 
LDQK  ITITLE) 
INIK1) 
INIK-1) 
ISKK12) 
LDO  ( IBLANK) 
ENIK14) 
LIUK  ISAVE). 
END 


ITLE»  LTI 
OF  6-WORD 
VERTED  TO 
NG  IN  FOR 
0202020B) 
ENIK1). 


LLS  (42). 
LLS  (6). 
SLJ  (1NEX). 
LLS  (42). 
STAKLTITLE) 


LLS  ( 
STAK 
LLS  ( 
SLJ  ( 
LLS  ( 
STAK 


42). 

LTITLE). 

6). 

2NEX). 

42). 

LTITLE). 


TLE) 
TITLE  RECORD. 
CODAP  IN  FORTRAN  62-3. 
TRAN  62-3. 
. 

SAVE  INDEX»  SET  COUNT. 

ENTER  02B. 

PERFORM 

SHIFTING. 
CHECK  IF  COMPLETE. 
COMPLETE  LAST  WORD. 
STORE  LAST. 
REPEAT 
FOR 

SECOND 
TITLE 
LINE. 


RESTORE  INDEX. 


FUNCTION  IPACK12  ( IONE ♦ 12 • I3» 14 ) 


PACKS  FOUR  12-BIT  WORDS 
WILL  REQUIRE  CONVERSION 


INTO  ONE  48-BIT  WORD. 
TO  CODAP  IN  FORTRAN  62-3 


LDA( I ONE) 

OLS(36) 

LDQ(I3) 

LLS(12) 

OLS(36) 

STAUPACK12) 

END 


LDO( 12) 
LLSU2) 
OLS(36) 
LDQ( 14) 
LLS(12) 


SUBROUTINE  IPACKL1  (LABELt  LABELK  LABEL2) 
PACKS  TWO  4-CHARACTER  LABELS. 
USE  DECODE/ENCODE  IN  FORTRAN  62-3. 
CONUFLAG  »  37773777B). 
LDA  (IFLAG)      LDQ  (LABEL). 
LLSC24)  STA  (LABEL1). 

LDA  (IFLAG)      LLS(24). 
STA  (LABEL2). 


00008450 
00008460 
00008470 
00008480 
00008490 
00008500 
00008510 
00008520 
00008530 
00008540 
00008550 
00008560 
00008570 
00008580 
00008590 
00008600 
00008610 
00008620 
00008630 
00008640 
00008650 
00008660 
00008670 
00008680 
00008690 
00008700 
00008710 
00008720 
00008730 
00008740 
00008750 
00008760 
00008770 
00008780 
00008790 
00008800 
00008810 
00008820 
00008830 
00008840 
00008850 
00008860 
00008870 
00008880 
00008890 
00008900 
00008910 
00008920 
00008930 
00008940 
00008950 
00008960 
00008970 


53 


END  00008980 

00008990 

FUNCTION  I CLOCK (I DUMMY)  00009000 

EXF  (01000B).  START  CLOCK.  00009010 

END  00009020 

END  00009030 


907 


■ 
54 


APPENDIX  III 
1.   Derivations  relating  ^   with  ^i  ,  and  $£     with  M     [VI. 

^x      2>X       <*K     W    L  J 

The  equation       £  =  ^7"  -f^^/^  /d£j ?) 
on  rearrangement  gives     .  "r     -  T^^A   /  — J  . 

re  srJ 

The  power  series   representation  of  the   inverse  hyperbolic  tangent   for 


values  of     /4HS\2  </  / 


is 


•Uc'(>£\  =  *£  +  4/2iSf+  d/&£f+  i(Z£±)\  .. 


9"  S 

so  that  to  a  good  degree  of  accuracy  it  can  be  written 


and 

rewriting 

C 

The  i 

8           di 

Since  the  depth  may  be  considered  as  being  a  function  of  C  d  =  F(C)|, 
and  C  =  G(X,Y),  then  by  the  chain  rule 

id     y  k  id     M  H 

dX        Dc    2K  %     dY  ~    be    »Y 


also  since  d  is  an  explicit  function  of  x  and  y 


a     &  ltd  ,  u\ 

dX*  M'    / DC2  2>c   '     ' 
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a     is    m    Mi 

iXrf         iKiY  I  DC*         DC  ' 


iY>         dY2   12>C>   *"  yc\     . 


2.   Derivations  of  the  partial  derivatives  of  the  depth  function  with 
respect  to  X  and  Y. 

The  derivatives  may  be  computed  directly  from  the  equation  for  the 
quadric  surface 

Jim  A,   +  A,X  -h  A3Y  +  A+KY  +    ASK2  +  ACY2 

so   that 

&..    A,   +  A<Y  y.  ZASK 

0X 

$    =    At    +A+X   +  ZAiY 


and 


ns  -  ZAS      i/  B    r  w 
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APPENDIX   IV 
Summary  of  refraction  and  shoaling  relationships   for 
intermediate  depths   [5J  [6  J. 
The  following  discussion  applies   to  waves  of   small   steepness  where 
the  deep-water  wave  height  divided  by  the  deep-water  wave   length   is   less 
than   .005    (H~/L     K  .005).      In  all   cases   the  subscript   zero   refers  to  deep- 
water  parameters . 

The  wave  velocity  depends  upon  wave   length  and  upon  the  depth  of 
water: 


2v  K   TC  J 


where  d  is  the  depth,  and  T  is  the  period  of  the  wave. 

Waves  of  a  certain  period  curve  as  they  approach  the  shore  from  deep 
water  until,  theoretically,  they  are  perpendicular  to  the  beach  in  zero 
depth  of  water.  For  any  change  in  depth,  Snell's  law  determines  the  cur- 
vature of  the  ray.   It  must  intersect  a  contour  at  an  angle  determined  by 
Snell's  law  for  the  successive  changes  in  depth-  The  tangent  to  the  wave 
ray  must  make  an  angle,  0(  ,  with  a  perpendicular  to  the  contour  at  the 
point  where  the  ray  intersects  the  contour.   The  ray  must  curve  with  the 
change  in  depth  so  that  Snell's  law  is  satisfied  at  a  discrete  set  of 
points  given  by  the  intersection  of  the  ray  with  a  set  of  contours.   As 
shown  in  Figure  IV,  the  wave  ray  crosses  the  contour  corresponding  to  the 
wave  speed  C  .   The  tangent  to  the  wave  ray  makes  an  angle  0(     with  a  line 
drawn  perpendicular  to  the  smoothed  contour.   Since  the  wave  ray  is  con- 
tinuously changing  direction,  it  must  make  a  new  angle,  £(„»  with  the 
perpendicular  to  the  contour  corresponding  to  the  wave  speed  C„,  when  it 
reaches  that  contour.   The  change  in  angle  is  A°{.      Then  at  the  two  con- 
tours corresponding  to  wave  speeds,  C-^  and  C2,  Snell's  law  holds  ^ince 
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the  wave  crests  intersect  the  contours  at  the  correct  angles.   The  impor- 
tant point  is  that  the  two  ray  tangents  are  connected  by  an  arc  of  a 
circle  which  determines  the  exact  path  of  the  wave  ray  from  point  A  to 
point  B.   The  iteration  procedure  described  in  the  text  is  used  to  arrive 
at  this  result. 

The  assumption  behind  the  wave  height  calculation  is  that  for  steady 
state  conditions  energy  does  not  flow  across  orthogonals  and  that  none  is 
destroyed  by  friction.   Therefore,  the  power  between  orthogonals  is  as- 
sumed to  remain  constant.   The  mean  wave  energy  per  unit  surface  area 
equals: 


£*  ietH* 


where  0   is  the  density  of  the  water  and  H  is  the  wave  height.   According 
to  wave  theory  only  a  fraction  of  the  wave  energy  is  carried  forward  with 
the  wave  form  at  the  speed  C.   Then  the  mean  power  transmitted  between 
orthogonals  equals: 

P  =  r\EC  d£ 

where  n  is  the  fraction  of  energy  carried  forward  and  dl  is  distance  be- 
tween orthogonals.   The  numerical  value  of  n  approaches  %  in  deep  water 
and  approaches  one  in  shallow  water.   By  equating  the  energy  in  the  deep 
water  to  that  in  the  shallow  water,  the  ratio  is  formed: 

jC  -  J  J      I       L 

4    2  n    cTcc    7!/%t 

where  the  terms   are  defined   above.     This   can  be  written: 


Mo  1  Z  /I     c/c.  < 
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where  K  is  the  refraction  coefficient  and  is  equal  to  -J  dl  /dl.   The 
term  under  the  square  root  sign  is  termed  the  shoaling  factor,  H  : 


//5  =/fiJZT 
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